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Symbols 


1 Tube flow 

2 Orifice flow 

A Area, ft 2 

a Angle-of-attack, deg 

C p Constant pressure specific heat, BTU/lbm - ° R 

Cq Discharge coefficient 

C v Constant volume specific heat, BTU/lbm - °R 

f Friction factor 

<p Roll angle 

g Gravitational Constant (32. 174 lbm-ft/lbf-sec 2 ) 

h Specific enthalpy, BTU/lbm 

K Thermal conductivity of Air (BTU/ft-sec °R) 

L Tube length, ft 

M Mach Number 

m Mass, lbm 

m Mass flow rate, lbm/sec 

P Pressure, lbf/ft 2 

R Tube radius, ft 

R e Reynolds Number 

p Density, lbm/ft 2 

T Temperature, °R 

t Time, sec 

U Internal energy, BTUAbm 

V Volume, ft 3 

V Velocity, ft/sec 


Subscripts 

a Mean value 

e Exit conditions 

i Inlet conditions 

m Manifold 
n Tube or orifice number 

oo Free Stream conditions 
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PREDICTION OF PRESSURE FLUCTUATION 
IN SOUNDING ROCKETS AND MANIFOLDED 
RECOVERY SYSTEMS 


SUMMARY 

The determination of altitude by means of barometric sensors is used in 
sounding rocket applications. Consequently, a method for predicting the per- 
formance of such sensing systems is needed. Herein a method is developed for 
predicting the pressure-time response of a volume subjected to subsonic air flow 
through from one to four passages. The pressure calculation is based on one- 
dimensional gas flow with friction. 

In addition, a computer program has been developed which solves the dif- 
ferential equations using a self-starting predictor-corrector integration tech- 
nique. The input data required are the pressure sensing system dimensions, 
pressure forcing function(s) at the inlet port(s), and a trajectory over the time 
of analysis (altitude-velocity-time), if the forcing function is trajectory depend- 
ent. The program then computes the pres sure -temperature history of the gas in 
the manifold over the time interval specified. 


INTRODUCTION 

This analysis, undertaken at the time of the development of the Aerobee 350 
Recovery System, has led to the development of a set of equations that describe 
the pressure and temperature histories within a manifold connected to a pressure 
source or sources. The source pressure histories are assumed to be known 
either as a function of Mach Number or time. In addition, a computer program 
which solves these equations is presented. 

This manifold system is used to initiate an operating sequence of a rocket 
vehicle or payload recovery system at a preset altitude. The chief difficulty 
which the analytical method has to deal with is predicting the sensitivity of a 
given manifold system to sensing the true altitude in terms of the pressures ex- 
isting at the external surface of the vehicle during flight. Further, the barostat 
connection to the pressure source is a variable. Therefore, in order to achieve 
the greatest generality, flow from the vehicle exterior to the manifold is assumed 
to pass through either tubes or orifices. Following is the development of the tube 
flow equations. 
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Assumptions 

The following assumptions are made in developing the equations used in this 
analysis: 

1. The pressure, density and temperature are distributed evenly and in- 
stantaneously throughout the manifold. 

2. The pressure, density and temperature at the port(s) are known for all 
times. 

3. The specific heats, C v and Cp are constant. 

4. The volume of the manifold is much greater than the volume of any tube 
leading into it. 

5. Continuum flow exists throughout the system. 

6. Entrance effects have a negligible effect upon the tube flow. 

7. An approximate equation for compressible adiabatic flow with friction 
can be used to calculate a mean value for velocity, given the mean 
density. 

8. Mass continuity is satisfied; i.e. , no mass addition in the manifold 
other than from the tubes. 

9. Steady flow exists over the integration interval. 

10. The behavior of air can be closely approximated by treating it as a 
perfect gas. 


DEVELOPMENT OF EQUATIONS 


The following theoretical development considers the case of tube flow from 
the external vehicle surface into the manifold and then discusses the variations 
in procedure required to account for orifice flow. A sketch of the flow case and 
general nomenclature is given in Figure 1. 

Where: T m = Temperature in the manifold 

T i = Temperature of the gas at the inlet port 

T e = Temperature of the gas at the outlet port 

t = Time - sec 
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ROCKET OUTER SHELL ADIABATIC WALL 



Figure 1 
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and 


dU 


m 


dt 


h = CpT 


Cv + T m m m ) 


h fi 


v 2 


flow 


ow 


2g 


= h = r T 

"source '-p 1 source 


dU 


m 


Now, solving for — — in Equation 1 


dt 

dU, 


dt 


m _ . / . . Vi 

- mj Ihj + 


V P 


2g J l he + 2g 


+ Q 


Substituting Equations 6 and 7 into Equation 8 

C v _ (rfliTj - rh e T e )Cp + Q 

Solving for T m and collecting terms: 


HQN 5 
EQN 6 

EQN 7 


EQN 8 


EQN 9 


4, Cp 2 n rh n T n + Q CpT m m m 

1 m 


C v m m 


EQN 10 


where rh m - 2 rh n , i.e. , there are no mass changes in the manifold other than 
those introduced by the flow. 

The mass flow rate in the tubes is determined as follows. Using an approx- 
imate equation for compressible adiabatic flow with friction. 


P - P = 
M r m 




v. 


EQN 11 


where p a and V a are mean values and i may be replaced by e where applicable. 

p a is determined by taking the average of the densities of air at the end 
of the tube being analyzed and the air in the manifold. V a is then found from: 


V 


a 



EQN 12 
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Knowing p a and assuming a value for V a , a Reynolds number and friction 
factor f is calculated for each tube and then a new V a is calculated from Equation 
12. V a (calculated) is compared to V a (assumed) and if the relative difference is 
greater than one percent, the iteration is continued until the process converges. 
Next, f is calculated from the following empirical equations: 


f = 0.0008 + 0. 05525/Re 0 ' 237 

Re > 1 00000 

EQN 13 

f = 0. 0791/Re 0 ' 25 

1 185 < Re < 100000 

EQN 14 

f = 16/Re 

Re < 1185 

EQN 15 


The p a and V a so found constitute a p a V a couple which is used to compute 
the mass flow rate through the tube being analyzed. 

m a = p a V a A EQN 16 

dp _ ^manifold _ ^ ^tubes EQN 17 

dt Voljnanjfoid V°lmanifold 

Equation 16 is used to compute the mass flow rate, rh n , for each tube. Since 
the dp/dt and dT/dt of the manifold are now known, we may integrate numerically 
Equations 10 and 17. The numerical integration technique is described in 
Appendix C. 


Orifice Flow 

For orifice flow there is one significant difference in the preceding calcula- 
tions: The mass flow rate through each orifice is calculated using 

ra a = C q A (^) K EQN 18 

where C q may be input as a variable and p is the density at the higher pressure. 

For both tube flow and orifice flow the limiting condition of choked flow is 
assumed to occur at M = 1 and this condition is applied by the computer program 
to limit the flow rate when necessary. 

The preceding analysis takes into account those variables which are felt to 
be significant for the problem analyzed. In most instances, the equations are 
programmed in an expanded form to facilitate checking and to provide a source 
of documentation to the user who may wish to identify the equations in their 
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programmed form. The program documentation will be found in Appendix A with 
a complete explanation of program variables and its use. 


Results 


Results from three tests are used to substantiate this analysis. First, a 
drop test of the Aerobee 350 recovery body at the El Centro Range in California, 
a DOD Parachute Test Facility. Next, data is used from the payload recovery 
of flight 17.05 GT-GG. Finally an experiment (Reference 10) which determines 
the pressure drop of a volume through various length and diameter tubes. The 
predictions of the manifold program will be compared to these tests and will be 
seen to match the test results closely. 

Actual pressure data from the Aerobee 350 drop test was used as a forcing 
function in the computer simulation. Drogue deployment was set for 20,000 feet, 
a manifold pressure of 972 psf. Actual deployment took place at 18, 700 feetwhere 
the manifold pressure reached 1000 psf. Measured deployment time was 33.35 
seconds after test initiation while the computer simulation predicted deployment 
at 33.66 seconds. In addition, prior to the drop test a pre-flight analysis of the 
maximum range of deployment altitudes was performed for several possible test 
configurations using wind tunnel test data (Reference 12). Results predicted that 
for a desired 20,000 foot recovery initiation, drogue deployment would occur be- 
tween 22, 500 and 15,000 feet for all configurations. Figures 2a and 2b show the 
forcing pressures and computed manifold pressure for the drop test. Thus, these 
results show the applicability and accuracy of this analysis. 

Pressure data from flight 17.05 GT-GG was also used as the forcing function 
in a computer program simulation. Figure 3a shows the forcing pressures and 
the computed manifold pressure while Figure 3b shows a comparison of the meas- 
ured and calculated manifold pressures. As may be seen from Figure 3b, the 
computed manifold pressure very closely follows the measured manifold pressure 
with a slight lag. This slight lag in the prediction may be due to any or all of the 
following: Non-uniform tube diameter, internal roughness of the tube, bends in 
the tube, or the possibility that the nominal manifold system used for computation 
is not the same as the actual manifold system flown. 

Another simulation of test data is presented in Figures 4 and 5 which show 
a comparison between program computations and some results obtained in Ref- 
erence (10) . The applied pressure differential and time shown in Figures 4 and 
5 is the time to equilibrium for each applied pressure. For the experiments 
performed in Reference (10), the experimental error is of the order of ±15 per- 
cent so that the program computation is considered to be a reasonably accurate 
simulation of the test. 
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These comparisons of results from the Aerobee 350 drop test and Flight 
17.05 GT-GG with the computations of the manifold computer program show very 
good agreement and tend to confirm the applicability of the analytical method. 
Similarly, the computed response times compare favorably with the results ob- 
tained in Reference (10) which further confirms the applicability of the equations 
for this specific application. Naturally, the equations and the manifold program 
will be subject to continued testing over diversified conditions in order to in- 
crease the confidence level as well as the proven regions of applicability. 


CONCLUSIONS 

This analysis, coupled with a computer program for solving the equations, 
can be used to predict the response time of a recovery system manifold 
connected to pressure sources on the surface of a re-entry body. In addition, 
the design of new manifolded recovery systems may be undertaken with confi- 
dence. The prediction of pressure variation in any ascending rocket vehicle is 
possible with the proper choice of a pressure forcing function and other system 
parameters. In general, as long as the assumptions of the analysis are met, 
the response of a volume subjected to a fluctuating pressure connected either by 
a tube(s) or orifice (s) may be analyzed. 
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PRESSURE (PSF) 



Figure 3b. Measured & Predicted Pressure for 17.05 Recovery 
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PRESSURE DIFFERENTIAL - mmHg 



RESPONSE TIME 


Figure 4. Pressure Response Through a Tube 




PRESSURE DIFFERENCE - mmHg 



1 10 100 1000 

RESPONSE TIME 

Figure 5. Pressure Response Through a Tube 
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APPENDIX A 


A1 The purpose of the manifold program is to solve the set of equations pre- 
sented in the body of this report. The solution to these equations is accomplished 
in a series of steps culminating in the numerical integration of Equations 10 and 
17 from the body of this report. 


_ 


Cp 2 n^n^n + Q CpT^rii 


m 


m 


C v m 


m 


_d_ 

dt 


m 


m 


' m 


= Z 


m r 


n V 


m 


EQN 10 
EQN 17 


The program is set up in steps so that various configuration and environ- 
mental parameters may be included in the solution of the equations, either to- 
gether or separately. This flexibility allows isolation of effects of the individual 
parameters on the pressure system being analyzed. In order to include this 
flexibility and yet make the input as simple as possible two concepts have been 
used with regard to the input. The first is the default option which means simply 
that the basic options needed to analyze the simplest case will be selected auto- 
matically in lieu of instructions to the contrary. The second concept is that the 
user need specify only those options which most suit the user's analytical model 
and data. The remainder of Appendix A will be used to describe the input vari- 
ables and their use. 


A2 The following are program variables/input names and are defined prior to 
describing the input options. 
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FLOW CHART OF MANIFOLD PROGRAM 


MAIN 

CALL INPUT SUBROUTINE 

TRANSFER DATA FROM INPUT 
ARRAY TO PROGRAM NAMES. 

INPUT PRESSURE COEFFICIENT 
DATA. 

CALL PRESSURE SUBROUTINE 
TO BEGIN CALCULATION. 



READIN 


INPUT SUBROUTINE 



INITIALIZE CONSTANTS 
INITIALIZE INTEGRATION ROUTINE 

START CALCULATION 

CALL INTEGRATION 
OUTPUT 

RETURN TO MAIN 
WHEN FINISHED WITH PROBLEM 



SINTRP 


QUADRATIC INTERPOLATION 
FOR MOST TABLES 


CSUBP 


LINEAR INTERPOLATION FOR 
FOUR WAY CSUBP TABLE 


INTEGT 


INTEGRATION ROUTINES 
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PROGRAM VARIABLE 
NAME AND 
INPUT NAME 


DEFINITION 


ALT 

CQ 

FRENO 

NCQ 

INDPVR 

PRINT FQ 

BDUMP 

TINIT 

TEND 

NOPTS 

ALPHA 

ALPHAA 

TALPH 

CASES 


Altitude Table — Ft. 

Discharge Coefficient table for orifice flow. 

The independent variable for the CQ table. The func- 
tion may be either Reynolds Number, the square root 
of Reynolds Number or pressure ratio where 
P/P 0 > 1. 

The n um ber of pairs of points in the CQ-FRENO table. 

Indicates the type of independent variable FREN0 = 

0 = (ReNo)*, - 1 = Re No, 1 = P/P 0 (P/P 0 > !)• 

The time interval between printed output steps — 
should always be used, otherwise every step will be 
output. 

BDUMP is input equal to 1 for all name load input to 
be printed, otherwise it is neglected. 

Initial time. 

Time at which program is to stop. 

Number of points in the Altitude -Velocity-Time Table, 
i.e. , number of time points. 

Input option to have the angle-of-attack vary with time. 
ALPHA is input only if this option is to be implemented, 
with ALPHA = Number of points in the a - T table 
(maximum of 10) . 

Name of the angle-of-attack array. 

Name of time array associated with ALPHAA. 

Number of cases to be run. 
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PROGRAM VARIABLE 
NAME AND 
INPUT NAME 


DEFINITION 


VR 

NEWFMT 

OMEGA 

PM 

RO 

TM 

PLTCON 

IPLOT 

PLOTFQ 

DELT 

DELMAX 

IERCRT 

ERCRIT(l) 

ERCRIT(2) 


Manifold Volume — input as inches 3 , converted inter- 
nally to feet 3 . 

Used to trigger NAME LOAD subroutine to read in user 
input case label. 

Roll rate in-C PS (optional input) . 

Initial Pressure in the manifold — PSF. 

Initial air density in the manifold lbm/ft 3 . 

Initial air temperature in the manifold ° R. 

Note: Only one of RO or TM need be input; the other is 
calculated in the program. 

= 1 Enables user to continue plotting from one case to 
the next, i.e. , could also be used when SAVE = 1 is 
used. 

The array in which the variables to be plotted are 
specified. 

Plot every PLOTFQ number of points (see A7 for 
explanation) . 

Initial time step, default = 0.0625. 

Maximum step size allowed by the user; if DELMAX 
is not specified there is no step size limiter. It is not 
normally necessary to limit step size but should be 
used to ensure that a minimum number of points are 
computed for output, i.e. , that the integration step is 
not so large that very little output is received. 

Error criteria for the numerical integrator. 

integration— preset. 
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PROGRAM VARIABLE 
NAME AND 
INPUT NAME 


DEFINITION 


JMAX 

N0NLIN 

ICP 

KT 

N0Q 

RT 

SAVE 


OPTSEL 

and 

CPCNTL 

IXMAX 

IYMAX 

IZMAX 


Number of equations to be solved— preset. 

Non-Linear Integration Option— preset. 

The number of different CP tables to be input, i.e. , 
succeeding cases may use either the same or a dif- 
ferent CP table. 

The number of tubes or orifices connected to the mani- 
fold, default = 1, Max No. - 4. 

Option selector for tube or orifice flow default = tube, 
N0Q = 2 for orifice. 

Recovery Temperature factor for aero heating in 
boundary layer. (Usually around 0.9). 

SAVE = 1 saves PM, TM and t from previous case to 
use in next case. This would be used when there is 
too much pressure forcing function data to input into 
the data array for a one case run. This would enable 
the user to just input more tables for each succeeding 
run. However, SAVE must be re -initialized for each 
case where it is desired . In addition, TEND must be 
respecified for each succeeding data set. 

OPTSEL is used as one part of a two part option selec- 
tor. OPTSEL selects the type of pressure function op- 
tion. CPCNTL is the second part of the two part option 
selector. CPCNTL selects Mach Number or time and 
symmetric or non-symmetric tables CP table definition. 

Number of Mach numbers or times to be input. 

Number of angle s-of-attack to be input. 

Number of 0-C p pairs to be input. The C p table is the 
pressure forcing function and may be defined as any of 
the input options d through i. 
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PROGRAM VARIABLE 
NAME AND 
INPUT NAME 


DEFINITION 


PIPSIZ 


Dummy array containing pipe or orifice dimensions 
(inches) . 

PIPSIZ 1 = Radius (1) 

2 = Length (1) 

3 = Radius (2) 

4 = Length (2) 

5 = Radius (3) 

6 = Length (3) 

7 = Radius (4) 

8 = Length (4) 


where length = 0 for orifice flow 


A3 General Option Description 

The different options are the following: 

Type of System 
Tube flow 
Orifice flow 

Parameters that may be included in the analysis 
Trajectory input 
Temperature function at the port 
Type of Pressure forcing function 
Type of discharge coefficient for orifice flow 

Roll rate. Note: This option requires that the integration step be 

smaller than other options require and therefore takes more com- 
puter time. 

Output all nameloaded input data 
Plot the results 

Include angle -of-attack if required by pressure forcing function tables. 

The preceding options are possible functions that the user may include in 
his analysis. It is up to the user to define his problem within the program's 
framework. 
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Detailed Option List: 


OPTION 

NUMBER/ DESCRIPTION 

LETTER 

0 Tube flow. 

2 Orifice flow. 

a Trajectory input required, altitude, velocity, time, 

b Ambient temperature used at port. 

c Recovery temperature used at port. TR = TA (1+2*RT*M) . 


d 

e 

f 

g 


Pressure forcing function AP 

C p = = f (M, a, 0). 

H q 

Pressure forcing function < 

^ 1 n /l K /TV\ I 


Cp = -=f(M,a,0)^-> lj. 
Pressure forcing function C p = P t = f (M, a, 0), lb/ft 2 . 


Pressure forcing function c _ _AP _ f ( t a 0 ) 


h Pressure forcing function P t ^ / p i , 

C p = p- = f(t,a,0),^->l 

i Pressure forcing function C p = Pi = f (t, a, 0) (lb/ft 2 ), 

j Variable Discharge Coefficient CQ = f(Re ,/z ). 

k Variable Discharge Coefficient CQ = f(Re). 


1 

n 

o 


Variable Discharge Coefficient 

CQ 




Input constant roll rate for vehicle simulation, cps. 


Output all input data except C p table which is always output 
automatically. 
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OPTION 

NUMBER/ 

LETTER 


DESCRIPTION 


p Plot up to 10 dependent variables as a function of any other 

variable. 

q Input angle -of-attack as a function of time to be used by the C p 

table. This table is needed only when C p is also a function of 
angle -of -attack. 

r Input non-symmetric C p tables , i.e., 0 < 0 < 360. 

s Use a simple one way table input into the altitude time array (when 

no trajectory is needed) to solve for the pressure fluctuation in a 
manifold. For example, this could be used to analyze the response 
of a pressure measuring system used in an experiment or to deter- 
mine the frequency response in a transducer-tube measuring 
system. 

t Hold the manifold temperature constant to affect an isothermal 

solution to a problem. This is accomplished by inputting the de- 
sired temperature with a minus sign. 

u Save the manifold pressure and temperature and time to use in the 

next case, i. e. , internal initialization of next case when multiple 
C p tables are required because there is too much C p data to input 
to one case. 


A4 Selection of options for input, i. e. , the input name and value input to imple- 
ment each option. Default options are those options the program will select auto- 
matically if no other option selection is made. 


OPTION 

INPUT NAME(S) 

OPTION SELECTION OR INPUT VALUE 

0 

N0Q 

0 (default) 

2 

N0Q 

2 

a 

ALT, VEL, TF 

Input data 

b 

RT 

0 (default) 



OPTION 

INPUT NAME(S) 

OPTION SELECTION OR INPUT VALUE 

c 

RT 

Equal to the recovery temperature factor 

d 

OPTSEL, CPCNTL 

0, 0 (default) 

e 

OPTSEL, CPCNTL 

-1, o 

f 

OPTSEL, CPCNTL 

-2, 0 

g 

OPTSEL, CPCNTL 

0, -1 

h 

OPTSEL, CPCNTL 

-1, -1 

i 

OPTSEL, CPCNTL 

-2, -1 

j 

INDPVR 

0 (default) 

k 

INDPVR 

-1 

1 

INDPVR 

1 

n 

OMEGA 

Roll rate cps 

o 

BDUMP 

1 

P 

IPLOT, NCURVS, PLTCON 

See A7 

q 

ALPHA 

Number of pairs of points in table. Then 


the table is input under arrays. 


Note: If an option is a default selection it need not (but may) be included in the 
input stream. 


OPTION 

INPUT NAME (S) 

OPTION SELECTION OR INPUT VALUE 

r and d 

OPTSEL, CPCNTL 

0 , 1 

r and e 

OPTSEL, CPCNTL 

-1, 1 

r and f 

OPTSEL, CPCNTL 

-2, 1 

r and g 

OPTSEL, CPCNTL 

0, 2 
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OPTION 

INPUT NAME (S) 

OPTION SELECTION OR INPUT VALUE 

r and h 

OPTSEL, CPCNTL 

-1, 2 

r and i 

OPTSEL, CPCNTL 

-2, 2 

s 

OPTSEL, CPCNTL 

CO 

1 

o 

t 

TM 

input -TM 

u 

SAVE 

1 


The input names OPTSEL and CPCNTL are used jointly to define the type of 
pressure function being input to the program. For some options, i.e. , d, e, f, 
it would not be necessary to input CPCNTL since 0 is the default value for 


CPCNTL but it might also be input solely for the sake of clarifying the input. In 
addition, the type of pressure function may change from case to case if so spec- 
ified. Otherwise, it will remain the same as the previous case, i.e. , symmet- 
ric option for case 1, non- symmetric for case 2, etc. 


A5 Table Definitions 

Input Names will be set off by quotes. Each type of table which may be input 
will be discussed. Quadratic interpolation is used unless otherwise specified. 

No extrapolation of any tables is performed. The end points are used as the ex- 
treme values. All tables, except the C table, are input via the name load sub- 
routine, which is input using standard Fortran formats. The method of using the 
name load option is described in A6. 

Option a calls for "Altitude", "Velocity", "Time" tables representing a tra- 
jectory experienced by the vehicle being analyzed with altitude and velocity being 
taken at the same time. 

"ALT" is the name of the altitude table 

"VEL" is the name of the velocity table 

"TF" is the name of the time table 

Each table may contain a maximum of 50 points. "N0PTS" is input as the 
number of sets of ALT-VEL-TIME points. The tables used for options d-i are 
linearly interpolated. This pressure forcing function table is input under stand- 
ard FORTRAN input rules, i.e. , no name load is used. The reason is that it 
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would be inefficient to input the amount of data required for the C p table under a 
name load option. The tables are specified as follows : 

"IX" = Number of Mach Numbers or time 

"IY" = Number of angle s-of-attack per Mach Number or time 

"IZ" = Number of roll angle-C p pairs per angle-of-attack 

The input format for IX, IY and IZ is (315) on the first card of the C p table 
data set. IX, IY and IZ are fixed point, i.e. , no decimal point, and are right 
justified. The C p 's for any option are then input as follows (Format = 2E15. 8). 

Mach Number or time whichever is the independent variable 
Angle -of-attack 

Roll angle-C p as determined by user input selection. These are input as 
follows. Suppose IX = 2 IY = 2 IZ = 2 

M[ or tj 

a i 

01 > c p 
02* C p 

<h 

0.. Cp 
0 2 , C p 
M 2 or t 2 

a i 

0i » Cp 

02. Cp 
02 

01 » C p 

02 » Cp 

Etc. 

The sample problems show the placement of the data on cards and position 
in the input stream. The C p tables are input after all the name load data for a 
particular case has been given. Up to 10 Mach Numbers or times may be input. 
Up to 15 angles -of-attack (per independent variable) may be input. Up to 20 0- 
C p pairs may be input per angle-of-attack. If the user requires more space for 
Cp data then option u would be selected allowing the continuation of the C p table. 
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Options j-i, the variable discharge coefficient table is input using the same 
load subroutines. 

"FRENO" is the name of the independent variable 

"CQ" is the name of the dependent variable 

Where FRENO may be Reynolds Number, the square root of Reynolds Num- 
ber, or a pressure ratio greater than 1. CQ is a value between 0 and 1. Each 
variable is input separately under the name load subroutine. There must be the 
same number of points in each table. "NCQ" is the variable that specifies this 
option in the input stream and is also the name used to specify the number of 
pairs of points in this table. 

Option Q is the angle -of-attack variation table. 

The name of the independent variable is "TALPH" 

The name of the dependent variable is "ALPHAA" 

Where TALPH is flight time and ALPHAA is angle-of-attack in degrees. 
This option is implemented by inputting the variable name "ALPHA" with the 
number of pairs of points in this table. A maximum of ten points may be input 
into each array. 

Examples of table input will be given in the section on Input. 

Note: All quadratically interpolated tables must have at least three points 
(per variable) input. The linearly interpolated tables must have at 
least two points (per variable) input. 


A6 Input Description 

The data are input to the manifold program via a "name load" read-in sub- 
routine, i.e. , each piece of data has a name and is identified in the input stream 
by this name. This also applies to most array variables in which case the input 
name is just the array name without subscript. A complete description of such 
a subroutine may be found in Reference (8) . This section will describe the 
method of using this name loader and will include sample inputs. 
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Name Load Input Options : 


NAME LOAD 
OPTION NUMBER 


OPTION 


0 

1 

2 

3 

6 

7 


Read in variable names, and data. Return to calling pro- 
gram, i.e. , manifold main program. 

Option 0 plus read another option card. 

Read in an array name, and data. Return to calling 
program. 

Option 2 plus read another option card. 

Read in an array name specific array locations, and data. 
Return to calling program. 

Option 6 and read another option card. 


5 Any of the above option numbers prefixed by a 5, i. e. , 

50, 52, 56, etc., specifies that a format will be input for 
that option. The format for the variable name read-in and 
array name read-in must be specified separately. Also, 
unless respecified the last format input will apply for either 
succeeding variables or arrays. 


Note 1: IF NO FORMAT is specified initially, a standard format of 10G8.4 will 
be used by both name load and array load options. This means that to 
change the format, an option number must be prefixed by a 5. However, 
the standard format should be sufficient for most users. 


Note 2: All data used in the manifold program, when it is input via name load 

routine, is input in floating point. The only fixed point data input to the 
manifold program are the table sizes for the M, a, 0 or t , a, 0 tables 
which are loaded in regular FORTRAN formats as described in A5. 
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CARD INPUT DETAILS 


Data Item 
Option Number 


Number of Variable 
names (maximum of 
10 per input card) 

1 

6-10 


NEWFMT (=0 or blank 
for no label card, 

=2 to input a label) 

1 

11-15 

315 

Card containing 
Hollerith format 

2 

1-80 

If NEWFMT set 
equal to 2 

Variable names left 
justified 

4 

1-8, 9-16, 
etc. 

10A8 

Format for Variable 
name input data 

5 

1-80 

Any valid Fortran 
format, default = 
10G8.4 

Array name left ' 
justified 

7 

1-8 


Location at which 
to start loading 
array (necessary only 
if loading does not 
start in location 1) 

7 

9-15 

A8, 17 

Array Name Format 

8 


Any valid Fortran 


format, default = 
10G8.4 


Card 

Type 


Column(s) 


1-5 


Format (Per Card) 


Fixed Point on C ard 
Type 1, right justified 
in field 
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Card Input Details 

A standard format for both variable name data and array name data is pre- 
set in the Read-in Subroutine. These formats may be used for inputting data 
without specifying a format in the input stream, i.e. , leave out card types 5 and 
8 if the preset format is large enough for the data to be input. These formats 
should be large enough for most applications and will save the user time in set- 
ting up input data. The preset format is the same for both variable name data 
and array name data and is (10G8.4) which permits 10 items to a card in fields 
of 8 columns each. Appendix B has a sample case with a listing of the input data. 

The Example inputs illustrate how the various kinds of data are input and 
are in no particular sequence. Line 1 indicates option number, column 5, num- 
ber of variable names, columns 9 and 10, and the number 2 in column 15 used to 
indicate that a Hollerith label will be input. Line 2 is the Hollerith label to be 
read in for this case. On Line 3 are the variable names to be read in while Line 
4 contains the data associated with each name. Since the name format and data 
format are each 8 columns long, the data appears under its name, however this 
need not be the case if the data input format is altered. Lines 5, 6 and 7 show 
the use of an array name, in this case the plot array. A part of the pressure 
forcing function input is shown in Lines 8 through 18. Line 8 gives the number 
of times (or Mach Numbers depending on the option), angle s-of-attack and roll 
angle points to be input. Line 9 contains the first time point. Line 10 contains 
the first alpha. Lines 11, 12, and 13 contain the roll angle -pres sure points. Line 
14 has the second angle -of-attack. Lines 15, 16, and 17 contain the set of roll 
angle -pres sure points associated with that alpha. Line 18 begins the sequence 
again with the second time point. Note that this C p table is a non-symmetric or 
360 degree table and as such would be so specified in the input stream even though 
it is not shown in this particular example which is abstracted from the total input 
shown on a later listing. 


A7 IPLOT — Name of the input array in which the plotted values are specified. 
The following values may be plotted: 

Index N^ 016 

1 Time 

2 Manifold Pressure 

3 Port pressure 1 

4 Port pressure 2 

5 Port pressure 3 

6 Port pressure 4 
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GODDARD SPACE FLIGHT CENTER 

FORTRAN CODING RECORD 


program MANIFOLD PROGRAM 

PUNCHING 

INSTRUCTIONS 

GRAPHIC 








PAGE OF 

programmer JOHN LAUDADIO | date 

PUNCH 








CARD ELECTRO NUMBER* 



GSFC 8-1(12/ 69) 








Index 


Name 


7 

8 
9 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

To set up this plot option use the following procedure: 

Set: NCURVES = Number of curves (maximum of 10) . 

PLTCON = 1 to continue plot from one case to the next, i.e. , con- 
tinuation of related cases. However, if unrelated cases 
are run set PLTCON = 0. 

IPLOT (1) = independent variable. 

(2) = dependent variable No. 1 by number. 

(3) = dependent variable No. 2 by number. 

Etc. up to 10 dependent variables. 

Note: PLOTFQ: If PLOTFQ is equal to PRNTFO, every printed point will also 

be plotted. However, by specifying PLOTFQ equal to an in- 
teger number n, only every nth point will be plotted. Notice 
that PRNTFQ may be any non-integer but only if PLOTFQ = 
PRNTFQ will every printed point be plotted. 


UNITS 

Altitude of vehicle— ft 
Velocity of vehicle — ft/ sec 


A8 Program Output 

OUTPUT NAME 

Altitude 

Velocity 


Atmospheric pressure 

Altitude 

Velocity 

Dynamic pressure 
Mach Number 
Angle-of-attack 
Manifold temperature 
Reynolds Number through port 1 
Reynolds Number through port 2 
Reynolds Number through port 3 
Reynolds Number through port 4 
Discharge coefficient for port 1 
Discharge coefficient for port 2 
Discharge coefficient for port 3 
Discharge coefficient for port 4 
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OUTPUT NAME 


UNITS 


Time 

Manifold Press 
Pressure at each port 
Atmospheric Pressure 
Q0 

Mach Number 
Alpha 

Manifold Temperature 
Pressure coefficient at each port 
Roll rate 

Mass flow rate (through each port) 
Reynolds Number (through each port) 
Internal and external ss 
Temperature (at each port) 

Phase angle at each port 
Total mass change in manifold 

Velocity in each tube 
Friction factor in each tube 


Time 

lbf/ft 2 

lbf/ft 2 

lbf/ft 2 

Dynamic pressure— lbf/ft 2 

Vehicle velocity/speed of sound at 
ALT 

Angle -of-attack— degrees 
Degrees R (Rankine) 

Same as input 

CPS 

lbm/sec 

Based on tube diameter 

Speed of sound— ft/sec 

Degrees R (T 00 or T r ) depending on 
input, i. e. , if recovery factor is used 

In relation to free stream, degrees 

lbm from t 0 to t present (TINIT to 
TEND) 

Average gas velocity— ft/sec 
Average friction factor in tube 



APPENDIX B 
SAMPLE PROBLEM 

As a means of illustrating both the type of problem that can be solved and 
the input details, the following sample problem is presented. 

Example 1 input — 

Data needed for program: 

Port sizes and tube lengths 
Type of pressure forcing function 

Port locations (fore and aft) needed only if using computed pressures 
at the ports or abstracting data from wind tunnel reports 

Trajectory data 


Option Selection — 

Option Number Description 


0 

a 

i 

b 


o 

r 

(non-standard) 

P 


Tube flow 
Trajectory input 
Pressure = f (t, a, 0 ) 

Atmospheric ambient pressure at each port 
Output all input data 

Input 360 degrees C p Tables 

Plot option for Port pressure 1, 2, 3, plus manifold and 
ambient pressure versus time 


u 

(non-standard) 


Save PM, TM and t to use in the next case 


B-l 



Option Implementation 

Option Input Names (Variables) 

General Input VR, PM, TM, TINIT, TEND, PRNTFQ, KT, CASES 

1 PIPSIZ 

a ALT, VEL, TF, NOPTS 

b Default, no input needed 

i, r OPTSEL, CPCNTL, C p Table indices and table values, 

number of C p tables, ICP 

o BDUMP 

p IPLOT, NCURVS, MODE, IPSKIP, PLTCON 

u SAVE 

The following data sheet is a listing of the program input for Example 1 . 
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COLUMN 


0 10 20 30 


1 10 ' 21 ! 
MANIFOLD ANALYSIS' OF 


[VR 

10 . 


PM, 

1080. 

6| 

MODE 

C.j 

2 


INCURVS 

is. 

! 1 
OPTSEL 
- 2 . 

3 

I PLOT 

1. 

3 

•ALT ■ 

1 7008.9 

16365.9 16302 
15754'. 7 

3 


TM ! 
520. 


I PSK'I P 

5. 


CPCNTL 

2-1 


2 .. 

2 Li 

j ! 

16|94 3. 7: 

9 


VEL 
3 38.5' 
330.2: 
306.5 
3 

TF ; 
360. 
362. 
364. 

2 

PIPSIZ 

.0625 

10 

360.2 

0. 

0. 

120 . 

240. 

5. 

0. 

120 . 

240. 

360.4 


21 

33 

32 

2 l 


3. 


1687is .9 
I624j0 .2 


|8.5 

18.5 


3610.2 

3^2.2 


lOj. 

2 


3 38. 
326. 


360. 

362. 


.062 


1048.6 
1062.9 

1183.7 

104e.6 

1062.9 

1133.7 


0. 

0. 

120 . 

240. 

5. 

0. 

120 . 

240. 

360.6 

0. 

0. 

120 . 

240. 

5. 

0. 

120 . 

240. 

360.8 

0. 


1013.9 
1041 .6 
1192.3 

1013.9 

1041.6 

1192.3 


1040.6 
1062.9 

1141.2 

1048.6 
1062.9 
1141.2 


A EPOPEE 

t i n I r 

360.2 

S4VF 

1. 


4. 


1,6014.8. 
1:617 7 .(7 


3 37.9 
3j24 . 6 


360.6 
3j62. 6 


TO. 


40 


50 


FLT 17. 0| 
TEND 
36 2. 

ICP 

2 . 


16749.9 
16 115. ) 


3 3 7|. 1 
322U4 


36 01.8 
3621. 8 


.0425 

I 


t 5 /;? 5 / 7 Oj 1 I 
|RLH;MP PfijNT F 0 

1 . .01 

PLICIS3 i 

i. ! 


7. 


60 


C AS ES 

2 . 


70 


80 


KT 

3. 


16685.4 16)621.2! 166517.4 
16054. 7 H993.4I 15933.2 


3 36. 3 
329. 2 


341 . 
3&i3 . 


335.5 j 334 
31 7.8 I 315 


1.2 ! 361 


3613 . 2 


1 0, 


36 3 J 


1649? 

15373 


3 3 3.2 
3|l 2.3 


361.6 
16 3.6 


NuPiT S 

21 J 


9 16429.4 


2 15 81 


331| 
30 9 


13.4 


.9 

.6 


36 1.8 
36 3.8 
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0 . 

120 . 

240 . 

5 . 

0 . 

120 . 

240 . 

361 . 

0 . 

0 . 

120 . 

240 . 

5 . 

0 . 

120 . 

240 . 

361.2 

0 . 

0 . 

120 . 

240 . 

5 . 

0 . 

120 . 

240 . 

361.4 

0 . 

0 . 

120 . 

240 . 

5 . 

0 . 

120 . 

240 . 

361.6 

0 . 

0 . 

120 . 

240 . 

5 . 

0 . 

120 . 

240 . 

361.0 

0 . 

0 . 

120 . 

240 . 

5 . 

0 . 

120 . 

240 . 

362 . 

0 . 

0 . 

120 . 

240 . 

5 . 

0 . 

120 . 

240 . 


1109.2 

1071.4 

1111.4 

1 109.2 

1071.4 

1111.4 


1079 . 

1084.2 

1090.2 

1079 . 

1084.2 

1090.2 


1078.8 
11 18.3 

1098.7 

1078.8 
1118.3 
1098.7 


1048.6 

1165.2 

1111.4 

1048.6 

1165.2 

1111.4 


983 . 5 

1062.9 

962.8 

983.5 

1062.9 

962.8 


1057.2 

1097 . 

1031.5 

1057.2 

1097 . 

1031.5 


1078.8 

1032.9 
1179.4 

1078.8 

1032.9 
1179.4 
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0 

TEND 

364. 

1 


10 

362.2 

0 . 

2 

3 

0 . 


1087.6 

120 . 


1075.7 

240 . 

5. 


1226.2 

0 . 


1087.6 

120 . 


1075.7 

240 . 

362.4 

0 . 


1226.2 

0 . 


1074.5 

120 . 


1071.4 

240 . 

5. 


1204.8 

0 . 


1074.8 

120 . 


1071.4 

240 . 

362.6 

0 . 


1204.8 

0 . 


110 M . 2 

120 . 


1105.5 

240 . 

5. 


1226.2 

0 . 


1109.2 

120 . 


110 5.8 

240 . 

362.8 

0 . 


1226.2 

0 . 


1074.5 

120 . 


1097 . 

240 . 

5. 


1200. 7 

0 . 


1074.5 

120 . 


1097 . 

240 . 

363 . 

0 . 


1200.7 

0 . 


1109.2 

120 . 


1105.5 

240. 

5. 


1204.8 

0 . 


1109.2 

120 . 


1105.5 

240. 

363.2 

0 . 


1204.8 

0 . 


1078.8 

120 . 


1050 . 

240. 

5. 


1192.2 

0 . 


107 e .8 

120 . 


1050 . 

240. 

363.4 

0 . 


1192.2 

0 . 


1117.9 



120 . 

589,5 

240 . 

5 . 

1234.7 

0 . 

1117.9 

120 . 

599 . 5 

240 . 

363.6 

0 . 

1234.7 

0 . 

1117.9 

120 . 

1126.8 

240 . 

6 . 

1213.3 

0 . 

1117.9 

120 . 

1 126.8 

240 . 

363.8 

0 . 

1213.3 

0 . 

1109.2 

120 . 

1126.8 

240 . 

5 . 

1204.8 

C . 

1109.2 

120 . 

1126.8 

240 . 

364 . 

0 . 

1204.8 

0 . 

1100.6 

120 . 

1118.3 

240 . 

5 . 

1204.8 

0 . 

1 100.6 

120 . 

1118.3 

240 . 

1204.8 



APPENDIX C 


The numerical integration scheme used to solve the equations in this report 
is a sixth order predictor corrector with a Runga-Kutta starter. Step size is 
automatically computed and altered based on both the stability and truncation er- 
ror. This particular predictor-corrector set was chosen rather arbitrarily dur- 
ing program design and is still used because it has never failed to work properly 
for the manifold program. The particular numerical integration program used 
was designed to handle up to fifteen differential equations simultaneously with 
added facilities for changing predictor-corrector pairs and error criteria easily 
as well as allowing utility in inputting the differential equations. 

Documentation of the above mentioned numerical integration scheme is now 
in progress. 

The predictor -corrector pair used is the following: 

Yn + 1 =.P = Yn + ^ [55 Yn - 59Y(n - 1) + 37Y(n - 2) 

- 9 Y(n - 3)] + At 5 Y 1V 

At • • 

Yn + 1 = C = Yn + ^ [9 Y(n + 1) + 19 Y(n) - 5 Y(n - 1) 

+ Y(n - 2)] - yjjjAt 5 Y IV 
where Y IV is the fifth derivative of Y. 

The four point Runga-Kutta formula used for starting the integration follows: 

Kj = At x f (x, y) 

a c ( At K i \ 

K 2 = At x f \x + — , y + — J 

K 3 - At x f (x + — , y + — J 

K 4 = t f (x + At, y + K 3 ) 
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Y =-^(K! + 2K 2 + 2K 3 K 4 ) 


¥;;=¥!+ AY 
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APPENDIX D 
PROGRAM LISTING 


IP-/* A 



//NQJFLMFD JOB ( NQOOl 1 879K , C» A00027 , HOOOOl ) , 122 ,MSGLE VEL = ( 1 > 1 ) 

7/ EXEC PGM=IEFBR14,REGI0N=6K 

//SCRATCH DO UNI T=D I SK , V0L=SER=M2SCR6 , DSN=NCJFLMFD, DISP= (OLO, DELETE ) 

// EXEC FORTRANG 
//SOURCE. SYSIN DO * 

C JOHN F LAUDADIO 3/18/70 

REAL*8DICT< 300) 

C Z ILCH1 IS COMMON FOR THE MANIFCLO SUBROUTINES 

COMMON/ZILCHl/YXl 15 ) , YPR I ME ( 1 5 ) , E RCR I T I 2 > , P ( A ) , R AD ( A ) , AL ( 4 » . 

1ALT< 50) ,VEL< 50) , TF < 50 ) , CQ I 10 > ,FRENO( 10) ,CMACH ( 10) , ALF A( 10, 15) , 
2ALPHAAI 10) »TALPH< 10 ) , PRESSF I 6 ) , TCON , TTUBE , TH ICK , CPTUBE , ROTUBE, 

3FI ( 10,15,20) ,CPP I 10,15,20) ,T, DEL T.DELMAX, PM, VR, THETA, OMEGA, 
4THRC0N, PLTCON , F A Z ANG , ALPH A , CP , A A, TM, PRNTFO, T END, PHI , RT, INDP VR, 
5JMAX,N0NLIN,IERCR T,I XMAX.I YMAX ,1 Z MAX , KT , NOQ, NCPT S , NCQ, N ALPHA, 
6NCURVS ,PLOTFQ, I PLOT ! 15 ) , APLOT ( 15 ) 

C THE • B * ARRAY CONTAINS ALL INPUT QUANTITIES FROM READIN 

DIMENSION BI 300) 

C. THE NAME BLANKS IS A FILLER IN VACENT LOCATIONS THAT 

C MAY BE USED FOR NEW VARIABLE NAMES AS NEEDED 

DATA DICT/50*'ALT*,10*'CQ',IG**FREN0*,*NCQ*,*PRNTFQ*,*BDUMP*, 
l'TINIT* ,' TEND * , ' NOP T S * , • ALPHA* ,* CASES* ,* VR* , 

2 'OMEGA* , 'PM' , • THRCON' ,' PL TCON* ,* BLANKS* ,* NOC* ,* RO* ,*TM* , • INDPVR* , 

3* BLANKS' , 'NCURVS* , • BLANKS* ,' PLCTFQ* ,7** BLANKS* ,4** PRESSF* , * OPTSEL* 
A, 'CPCNTL ' ,' TCON* , * T TUBE* ,* THICK* ,* CPTUBE* ,* ROTUBE* ,5** BLANKS*, 

5 * DELT* ,'PELMAX' ,2*'ERCRIT* ,« JMAX* , * NCNLIN* ,* IERCRT* ,* ICP*,'KT*, 

6 * SAVE* , 10*' ALPHA A* ,10*' TALPH* ,* RT* ,4*' BLANKS* , 

750*' VEL' ,50*' TF* ,8**PIPSIZ*, 15**1 PLOT*, 27** BLANKS*/ 

ND ICT=300 
DO 12 1=1,300 
12 B I I ) =0 . 

200 CALL READINIDICT, B.NDICT ) 

55551 CONTINUE 
NCQ=B (71) 

PRNTFQ=B ( 72 ) 

BDUMP = B ( 7 3 ) 

C T IS THE INITIALIZATION TIME FOR THE INTEGRATION 
C T=T l N I T = INITIAL TIME 

T=8 ( 74 ) 

TEND=B ( 75 ) 

C NOPTS = THE NUMBER OF ALT-VEL TF DATA POINTS INPUT TO THE PROG 

NOPTS=B ( 76 ) 

ALPHA=8 ( 77 ) 

ICASES=B(78) 

VR=B( 79) 

OMEGA=B ( 80 ) 

PM=B(81) 

C THERMAL CONSTANT K IN BTU/SEC*FT*DEG R 

THRCON=B« 82) 

C PL TCON SET =1 WILL ENABLE THE USER TO CONTINUE PLOTTING EROM 

C ONE CASE TO THE NEXT, WHEN THE LAST CASE TO BE PLOTTED IS 

C REACHED SET PL TC0N=0 , OTHERWl SE PLTCON IS SET = 0 AUTOMATICALLY 

C WHEN THE LAST CASE IS INPUT 

PLTC0N=8 ( 83 ) 

I F ( ICASES .LF. 1) PLTC0N=0. 

B ( 83 ) = PLTCON 

C OPTION SELECTION , 0=TU8E FLCW,l= TUBE FLOW+H EAT, 2= ORIFICE FLOW 

C IF NOQ IS NOT INPUT THE TUBE OPTION IS CHOSEN 

N0Q = 8 I 85 ) 

C YX ( l ) IS THE INITIAL DENSITY-DERIVED FROM INITIAL TMGPM IN 

C SUBROUTINE PRESUR UNLESS DENSITY IS INPUT ,IE DENSITY IS 
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c NONSTANOARD AND TAKES PRECIOENCE IF BOTH TEMP AND DENSITY ARE 

C INPUT BY MISTAKE. PRESS MUST ALWAYS BE INPUT 

C YXI2I IS THE INITIAL TEMPERATURE 

YX( I ) =B I 86) 

YX ( 2 ) =B ( 87 ) 

I F I YX( 2 ) .LT. 0.) WRITEI6, 10001 

1000 FORMAT! *0 THE ISOTHERMAL SOLUTION TO THE EQNS HAS BEEN CHOSEN*) 
INOPVR=B( 88) 

C INDPVR IS A TRIGGER INDICATING THE TYPE OF INDEPENDENT VARIABLE 

C BEING USED FOR THE CQ TABLE, IE, -1,REN0, OIDEFAULT ) SORT (RE) » 

C l PRESSURE RATIO .GT. .1.0 

IF ( I NDPVR .NE. 0) WRI TE ( 6 , 100 ) I NDPVR 
100 FORMAT! *0 A NONSTANDARD INDEPENDENT TABLE INPUT HAS BEEN USED FOR 
ICQ, INDPVR=',12) 

C SET UP VARIABLES FOR PLOTTING 

NPPTS=B ( 89 ) 

NCURVS=B ( 90 ) 

C DEFAULT NPPTS 

IF(NPPTS .EQ. 0 .AND. NCURVS .GT. C) NPPTS=100 
MC0E=B(91 ) 

PL0TFQ=8 ( 92 ) 

DO 205 l X= l » 1 5 
APLOT( IX)=B(258+IX) 

205 IPL0T(IXJ=B(258+IX) 

C NOTE THAT PRESSF(5C6) ARE INPUT AS VARIABLE NAMES IOPTSEL C CPCNTL) IN 
C READIN BUT IS HANDLED EVERYWHERE ELSE AS PRESSF (5&6 ) 

C THE NAME INPUT, CPCNTL, IS USED TO MAKE INPUT EASIER 
DO 206 I X= 1,6 

206 PRESSF! IX)=B(99+IX> 

I F (PRESSF ( 5 ) .EO. -1.) WRITE(6,110) 

110 FORMAT! *0 THE SPECIAL USE OPTION HAS BEEN CHOSEN, IE, P(K)=PO*CP*) 
IF ( PRE SSF ( 5 ) .EO. -2.) WRITE(6,ll3) 

113 FORMAT! *0 THE SPFCIAL USE OPTION HAS BEEN CHOSEN, IE, P(K)=CP ') 

IF ( PRESSF ( 5) .EO. -3.) WRI TE (6,111 ) 

111 FORMAT! *0- THE SPECIAL USE OPTION HAS BEEN CHOSEN, IE, P(K)=F(TF)*) 
IF ( PRESSF ( 6) .EO. -1. .OR. PRESSF (6 ) .EC. 2. ) WRI TE (6 , 112 ) 

112 FORMAT COTHE SPECIAL USE OPTION OF CP AS A FUNCTION OF TIME HAS 
1BEEN CHOSEN, IE, CP=F ( T , ALPHA ,PHI ) * ) 

C FOR THE SPECIAL USE OPTIONS THE ATMOSPHERIC VARIABLES ARE INPUT 

C AS CONSTANTS IN THE PRESSF ARRAY AS FOLLOWS 

C RHOO=PRESSF(l) , PO=PRESSF (2) ,T0=PRESSF(3), A8=PRESSF ( A )= SOUND SPD 

C THIS DEFINES THE ATMOSPHERIC QUANTITIES WHEN THEIR CALCULATION 

C IS PRECLUDED BY THE SPECIAL USE OPTIONS 

C PRESSF (6) =-l. ,0., FOR SYMMETRIC CP TABLES AS FIT) »F (MACH NO), 

C PRESSF (6) =1 ,2. , FOR NGN-SYMMETRIC CP T ABL ES,CP=F ( MACH ) ,CP=F ( T) , 

C SYMMETRIC TABLE = 180 DEG , NONSYMMETRIC =360 DEG INPUT TABLES 

TCCN =B( 106) 

TTU8E=B( 107) 

THICK=B( 108)/12. 

CPTUBE=8( 109) 

R0TU8E=B l 1 10) 

C B ARRAY LOCATIONS 111-115 INCLUSIVE ARE AVAILABLE FOR USE 

DELT=B( 116) 

C DEFAULT DELT IF IT IS NOT SPECIFIED 

I F ( DELT .EO. 0.) CELT=.0625 

C IF HEAT OPTION IS CHOSEN OELMAX WILL BE SET EQUAL TO .001 

C UNLESS OELMAX IS ALREADY SET LESS THAN .001, IE, 

I F ( NOQ .EQ. 1 .AND. (8(117) .GT. .001 .OR. B( 1 17 ) . EQ.O . ) )B( 117 )= .001 
DELMAX=B( 117) 

ERCRIT! 1 )=B< 118) 
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c 

c 


c 

c 

c 


c 

c 


105 

106 
C 


225 

C 

C 


PRESENT CASE TO BE USED IN NEXT CASE 
TIME IT IS TO BE USED 


226 


107 


10001 

C 

c 


£RCRIT(2)=B( 119) 

JMAX=B( 120) 

DEFAULT JMAX , NOTE JMAX IS NEVER INPUT UNLESS THE NO OF 
EQUATIONS IS CHANGED 
I F ( JMAX .EO. 0) JMAX=2 
IFIBI120) .NE. JMAX)B(120)=JMAX 
NCNL I N=0 (121) 

IERCRT=B( 122) 

ICP DENOTES THE NUMBER OF CP TABLES TO BE USED DURING A RUN 
ICP=B( 123) 

KT IS THE NUMBER OF PURTS ENTERING THE MANIFOLO 
KT=B( 124) 

DEFAULT KT= l 
IF (KT .EQ. 0) KT = 1 
I F ( B ( 124 ) .NE. KT) B(124)=KT 
SAVE=1 SAVES PM »PM , T FROM 
SAVE MUST BE RESET EACH 
SAVE=B( 125) 

I F ( NOQ .LE. 1) WRITE (6,105) NOQ 
I F ( NOQ .GE. 2) WRITE (6,106) NCQ 
FORMAT ( ' 0 TUBE FLOW HAS BEEN CHOSEN, N0Q=*,I3) 

FORMAT ( * 0 ORIFICE FLOW HAS BEEN CHCSEN, N00=',I3) 

LOAD ALT.VEL, TIME ARRAYS 
DO 1 1=1,50 
ALT ( I) =B ( I ) 

VEL ( I )=B( 1+150) 

TF( I )=B( I+2C0) 

I F ( ( ALT ( 1 ) .NE. 0. .OR. VEL(l) .NE. 0. .OR. TF(1) .NE. O.J.AND. 
1NCPTS .EQ. 0.) WRITE(6,225) 

FORMAT ( • 0 THF TRAJECTORY HAS BEEN INPUT BUT NOPTS IS STILL REQD* ) 
LOAD ORIFICE COEFF ARRAY AS A FUNCTION OF RENO, IE, FRENO 
LOAD ANGLE OF ATTACK -TIME TABLE 
DO 5 1=1,10 
CQ( I )=B( 1+50) 

FRENOI I )=B( 1+60) 

ALPHA A ( I ) =B( 1+125 ) 

TALPHI I ) =8 ( 1+135) _ „ , 

IF(ALPHAA(l) .NE. 0. .OR. TALPH(l) .NE. 0. .AND. ALPHA .EQ. 0.) 

1 WRITE (6, 226) 

FORMAT! *0 ALPHAA OR TALPH HAS BFEN INPUT BUT ALPHA IS STILL REQD') 
1F(C0(1) -NE. 0. .OR. FRENOI 1 ) .NE. 0. .AND. NCC .EC. 0.)WRITE(6, 

Format ( • o orifice coeff have been loaded but nco (count of ncq"S) 

1HAS NOT BEEN LOADED* ) 

I F ( BDUMP .EO. 1.) WRI TE (6 , 10001 ) ( B( I ) ,1=1,300) 

FORMAT ( IHO , 1 0G1 3. 6 ) 

RT GT 0 IMPLEMENTS THE RECOVERY TEMP OPTION , RT IS THE RECOVERY 
TEMPERATURE FACTOR .USUALLY ABOUT .9 


RT=B( 146) 

WRITE ( 6 , 906 ) 

906 FORMAT ( IHO, ' INPUT OIM, VR IN IN., OMEGA IN CPS, RAD(K) IN IN., 
1AL(K) IN IN., ALL OTHERS ARE IN FT, LBF , LBM,S EC' ) 

I F ( BOUMP .EQ. 0) WRITE(6, 902) VR, OMEGA, PM, YX(2),TINIT 

902 FORMAT( IHO, ' INI TI AL CONDITIONS ' , ' VR=* , G13 .6, • OMEGA= • , G1 3. 6, * PM= ' , 
l G13. 6, • TM=* ,G13.6,'TINIT=* ,G13.6) 

C CP TABLES INPUT HERE 

I F ( ICP .EQ. 0) GO TO 921 
ICP=ICP-l 
B(123)=ICP 

READ ( 5, 950) I XMAX , I YMAX , I ZMA X 


D-3 



O O >0 <-u o o o 


950 FGRMAT (415) 

DC 920 I X=1 * I XMAX 
RE AO ! 5, 9 l 5 ) CMACH(IX) 

WR l TE ! 6 , 904 ) CMACHUX) 

904 FORMAT! 1H0.8G10.5 ) 

915 FORMAT ( E15 . 8 ) 

DO 920 I Y= 1 * 1 YM AX 

READ! 5,916) ! ALF A ! I X ,1 Y ) , ( F l( IX , I Y , 1 1 ) ,CPP( IX, I Y, 1 Z > , I Z=l, IZMAX ) > 

WRITE! 6,916) < ALF A ( I X , I Y ) , ( F I ( I X , I Y , I Z ) ,CPPI I X , I Y, I Z ) , I Z=1 , I ZMAX ) ) 

916 F0RMAT!E15.8/!2E15.8) ) 

920 CONTINUE 

921 CONTINUE 

PIPE SIZE (PIPSIZ ARRAY) INPUT FROM 8(250-258) 

PIPSIZ DOFS NOT APPEAR AS AN ARRAY NAME IN THE PROGRAM 
IT IS ONLY A MEMNCNIC FOR DATA INPUT 

1 = 1 

DO 305 K=2 ,8,2 
RAO ( I )=B ( 249+K) 

AL ( I ) =B ( 250+K ) 

1 = 1 + 1 

05 CONTINUE 

VR=VR/1728. 

DO 92 2 11=1,4 
RAD! I I ) = R AO ( II )/12. 

22 AL! 1 1 ) = AL( 1 1 )/12. 

INITIALIZE PHI HERE IN CASE A NON STANDARD DPTION IS USED 
PHI IS USEO IN CALCULATING THE RCLL ANGLE IN CSUBP 
PH 1 = 0. 

14 CCNTINUE 

C INITIALIZE MANIFULO SUBROUTINE HERE 

C ENTER CALCULATION PHASE OF MANIFOLD PROGRAM 

CALL PRESUR 

I F ( SAVE .NE. 1.) GO TO 55552 
B ( 81 ) =PM 
B ( 8 7 ) = TM 
B ( 74 ) =T 
SAVE=0. 

55552 CCNTINUE 

ICASES=ICASE S-l 
B ( 78 ) =ICASES 

I F ( ICASES .GT. 0) GO TO 200 

IFdCASES ,LE. 0 .ANO. ICP.GF.i) GO TO 55551 

STOP 

END 
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SUBROUTINE PRESUR 

COMMON/ZILCHl/YXI 15) f YPR I HE ( 15 > , ERCRIT ( 2 ) , PI 4 ) , R AO ( 4 ) , AL 1 4 ) . 
1ALTI50) ,VEL(50> , T F I 50 ) , CQ 1 10 ) , FRENO 1 10 ) .CMACHI 10), ALFA! 10,15), 
2ALPHAA(10),TALPH( 10 ) , PRESSF ( 6) , TCON, TTUBE .THICK, CPTUBE, ROTUBE, 

3FII 10,15,20) ,CPP< 10,15,20) , T ,OELT .OELMAX , PM, VR, THETA, 0* EGA, 
4THRC0N* PLTCON , F A Z ANG , ALPHA, CP, A A,TM, P RNTFO, TEND, PHI, RT, INDPVR, 

5 JMAX , NONL IN, IERCRT , I XMAX, IYMAX, lZMAX,KT,NOQ,NOPTS,NCQ,NALPHA, 

6NCUR VS, PLOTFQ, I PLOT I 15) ,APL0T(15) 

INTEGER ONE/1/, TWO/2/, 1/58/ 

C OP/OT IN THIS PROG IS BASED ON ENERGY AND HEAT TRANS EQNS 

DIMENSION ANSI 4) ,0M14) ,TAI4) ,RH0I4) ,C(15) ,CCWI4) ,CPW 14 ) , REW ( 4 ) 
DIMENSION PHI I (4) ,PLOT < 22 ) , VVOUT 1 4) , QAI R 1 4 ) 

DIMENSION TUBE (4) , TCONV ( 4 ) , AVOL ( 4 ) , AL AT I 4 ) , TVOL ( 4 ) 

DATA ISKIP/O/ 

INITIALIZE CONSTANTS FOR CALCULATION 

DEFAULT FOR ORIFICE EQNS IF TUBE DIMENSIONS ARE ZERO 
I F ( AL ( 1 ) .EQ. 0. .AND. NOQ .LT. 2) N0Q=2 

C INITIALIZE V2 TO DETERMINE CPTICN SELECTION 

V2=0. 

C INITIALIZE PLOT ROUTINE HERE 

I F ( NCURVS .GT. 0 .AND. ISKIP .EQ. 0) CALL RJPLOT ( APLOT, NCURVS) 

C IF PL TCON IS USED FOR THIS CASE SKIP PLOT REINITIALIZING 

I F ( PLTCON .EQ. 1.) ISKIP=1 
I F ( PLTCON .EQ. 0.) ISKIP=0 
I IPPSS=0 
P I =3. 14159 
GAMMA= 1 . 4 
VV = 1 • 

A J=778. 26 

C OMEG=ROLL RATE IN CPS 

C SAVE OMFGA IN CPS THEN CONVERT TO RAD/SEC 

OMEG=OMEG A 
0MEGA=0MEGA*6. 28318 
GC=32. 174048 

C R IN ( FT-L8F ) / ( LBM- DEG R) 

R = 5 3.36 
CCP = • 240 
CV= .1710 

C INITIALIZE TOTMAS WHICH IS THE TOTAL MASS CHANGE IN THE MANIFOLD 

TOTMA S = 0 . 

RR = ( ( .75*VR)/PI )**( .3333) 

TM=ABS( YXI 2) ) 

C TISO IS THE TRIGGER FOR KEEPING THE TEMP CONST ,IE ISOTHERMAL 

T I SO=YX( 2 ) 

C MAKE SURE YX 1 2 ) IS POSITIVE FOR INTEGRATOR 

YXI2)=ABS(YX(2) ) 

RO=YX( 1) 

C IF YX I 1 ) NE 0 USE DENSITY INPUT OPTION 

IFIYXI 1 ) .EQ. 0. ) GO TO 70 
TM=PM/(R*RO) 

YX ( 2 ) =TM 
GO TO 75 

70 RO=PM/lR*TM) 

YXI 1 ) =RQ 

75 CONTINUE 

C ZERO OUT PORT VARIABLES 

00 11111 K= 1,4 
P ( K ) = 0. 

T A I K ) =0 . 

CPW I K ) =0 . 
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DM ( K ) =0 . 

REW ( K ) =0. 

PH 1 1 ( K ) =0. 

VVCUT ( K ) =0. 

CGW ( K ) =0. 

C INITIALIZE INITIAL HEAT TRANSFER RATE AND CALCULATE CONSTANTS 

QAIR(K)=0. 

TUBE I K ) =TTU8E 
TCCNV < K ) = TC0N 
AVOL<K)=PI*RAD(K) **2*AL(K> 

ALAT<K)=2.*PI*RAD(K>*AL(K) 

TVCL(K)=PI*( RADIK ) +TH ICK ) **2*AL < K J-AVOL IK) 

WRITE (6, 1010<J)AV0L(K) , AL AT < K > ,TVOL I K ) 

11111 RHO I K ) =0. 

C SET V2=PRE SSF I 5 ) IF PRESSFI5) LT 0 

I F ( PRESSF ( 5 ) .LT. 0.) V2=PRESSF(5> 

C INPUT ATMOSPHERIC VARIABLES HERE IF V 2 IS NEGATIVE 

IFIV2 .EO. 0.) GO TO 81 
RHOO= PRESSF 1 1 ) 

PG=PRESSF ( 2 ) 

TO=PRE SSF ( 3 ) 

AB=PRESSF I A ) 

AA=0. 

H0=0. 

V1=C. 

00 = 0 . 

CP=0. 

81 CCNTINUE 

C FAZANG=ANGLE between kt evenly spaced ports 

FAZANG=(2.*PI )/KT 

C DISCHARGE COEFFICIENT FOR AN ORIFICE -CONSTANT 

CGD=.611 

C INITIALIZE INTEGRATION SUBROUTINE HERE 

C LIMIT INTEGRATION STEP SIZE TO OUTPUT FREQUENCY 

I F ( DELMAX .GT. PRNTFQ ) DELMAX=PRNTFQ 

CALL START 1 (YX,YPRIME,T,UELT,DEIMAX, ERCRIT, JMAX.NONLIN, IERCRT) 
ALPHA INITIALLY IS THE NO OF POINTS IN THE ALPHAA-TALPH 
TABLE AND IS USED THEREAFTER TO TRANSFER ALPHA TO CSUBP 
I ALF IS THE NO OF POINTS IN ALPHAA TABLE 
I ALF=ALPHA 

THE DERIVATIVES ROD ANO TMD ARE CALCULATED IN PCALC 
IS=0 
1 1 1 = 1 
I0RDER=3 

PLOTF IS INITIALIZED HERE 
PLCTF=0. 

5 CCNTINUE 

TPRNT IS THE TIME AT THE LAST OUTPUT STEP AND IS SAVEO TO 
COMPARE WITH PRNTFQ I PRINT FREQUENCY) 

TPRNT =T 
l CCNTINUE 

FINO ALT C VEL AT TIME T 
I F ( I ALF . EO. 0 ) GO TO 45 

NOPRNT = 3 SUPPRESSES PRINTEO OUTPUT FROM SINTRP WHEN TABLE 
BOUNDARIES ARE EXCEEDED, SINTRP CHOOSES END POINTS 
N0PRNT=3 

CALL SINTRP! T,ALPHAA,TALPH»I ALF, I S»II I , I ORDER, ALPHA, NOPRNT ) 

45 CCNTINUE 

N0PRNT=3 

I F (NOPTS .EO. 0 > GO TO 80 
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CALL SINTRPI T.ALT.TF ,NOP T S , I S , l II , I CRDER,HO, NCPRNT > 

NGPRNT=3 

I F ( PRESSF ( 5 ) .EQ. -3.) GC TO 80 

C SKIP AT62 ANO VEL INTERPOLATION WHEN 2ND SPECIAL USE OPTION 

C IS USED 

CALL SINTRPIT,VEL,TF,NOPTS, IS, III , IORDER, Vl.NOPRNT > 

CALL AT62IHO ,ANS) 

RHOO=ANS( 1) ♦GC 
PO=ANS( 2 ) 

T0=ANS<3)*1. 80007 
AB=ANS( A) 

A A= ABS ( VI) / AB 
QO=(RHOO*V1*VI*.5)/GC 
80 CONTINUE 

30 CONTINUE 

C BEGIN CALCULATION OF DERIVATIVES RODIRO DOT) £ TMDITM DOT) 

THETA=C. 

DMT0T=0. 

SUM=0. 

SUMP=0. 

PM=R*RO*TM 

C INITIALIZE HEAT TRANS TERMS BEFORE LOOP 

QTCT=0. 

DTA=0. 

DTT=0. 

10003 DO 500 K= 1 » KT 

C V2=-l , P ( K ) = PO*CP , V2=— 2 , P I K ) =CP , V2=-3 ,PO=FCTF) WITH 

IC0NS=2 
COUNT=0. 

I F ( V2 .GE.-2.) CALL CSUBP 
P ( K ) = PO + C P*QO 

C I F ( PRESSF ( 5 ) LT 0) PORT PRESSURE IS ALTERED FOR SPECIAL USE OPTION 
C THE FORCING FUNCTION INPUT INTO THE ALTITUDE ARRAY 

IFIV2 .EQ. -1.) P(K)=PO*CP 
IFIV2 .EQ. -2.) P ( K ) = CP 
I F I V 2 .EQ.--3.) P ( K ) =H0 

C IE P(K) GOES TO 0. THE DENSITY ECN WILL BE DIVICING BY 0. ALSO 

IFIP(K) .LT. 0.) P ( K ) =0 . 

C RT IS THE RECOVERY TEMPERATURE FACTOR 

1020 T A ( K ) =T0* ( l.+.2*RT*AA*AA) 

RHO(K)=P(K)/(R*TA(K) ) 

C TE IS THE TEMP USED TO CALC THE ENERGY TERM OF DT/DTIK) 

C USE TEMP OF SOURCE GAS FOR TF 

TE = T AIK) 

IF(PIK) .LT. PM) TE=TK 

C USE AVG DENSITY FOR APPROXIMATE COMPRESSIBLE ADIABATIC 

C FLOW WITH FRICTION THROUGH A PIPE, IE ASSUMMED LINEAR VARIATION 

C USE DENSITY FROM DIRECTION OF FLOW FOR ORIFICE FLOW 

C CALCULATE THE SPEED OF SOUNO -SS IFT/SEC) BASED ON THE 

C TEMPERATURE FROM THE DIRECTION CF FLOW 

SS=SQRT(GAMMA*GC*R*TE) 

C CALCULATE VISCOSCITY BASED ON TE 

AKUU= 2.22E-8*(TE**.5/(1.*(180./TE))) 

C IF NOC= 2 OR 3 USE ORIFICE FLOW EQNS 2 NO HEAT , 

IFINOQ .LT. 2 .OR. NOQ .GT. 3) GO TO 601 
ROAVG=RHO!K) 

IF(PIK) .LT. PM) ROAVG=RO 
DMSAVE=0. 

C ROAVGS=DENSI TY IN SLUGS ,CQD = DISCHARGE COEFFICIENT 

C ROAVG =DENSl TY IN LBM/F T**3 OMK) IN LBM/SEC 
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RQAVGS=RO AVG/GC 

C ORIFICE FLOW EQUATIONS 

AO CONTINUE 

DM ( K ) =CQD*PI *RAD( K) **2*ROAVG*SCRT ( 2.*( ABS (P(K )— PM) ) /ROAVGS ) 
IF(PIK) .LT. PM) OM IK) = -DM ( K ) 

VTR I L=ABS ( VV ) 

VV=OM(K)/ ( ROAVG*P 1*RAD(K)**2) 

RE=(ROAVG*ABS(VV)*2.*RAO(K) )/(AMUU*GC> 

I F ( VV .GT. SS .AND. NCQ .EQ. 0) GO TO 600 
IF (NCQ .EO. 0) GO TO 1060 
I F ( VV .NE. 0.) GO TO 50 
VV=1. 

DM(K)=0. 

GO TO 1060 
50 CONTINUE 

C ITERATE TO FINO OM(K) IF COO IS VARIABLE 

VTEST=ABS(VV/VTR1L) 

IF(ABS(VTEST-1. ) .LE. .01 .AND. VV.GT. SS) GO TO 600 
IF< ABSI VTEST-1. ) .LE. .01) GO TO 1060 
C THE DISCHARGE COEFF TABLE IS CQD VS (RE)**. 5 

C CHECK TO FIND WHICH I NDE P VAR IS USED FOR CQ 

IF( INDPVR)76,77,78 

76 RS=RE 

GO TO 79 

78 RS-P ( K ) /PM 

I F ( PM .GT.P(K)) RS=PM/PIK) 

GO TO 79 

77 RS=SQRT ( RE ) 

79 CONTINUE 
N0PRNT=3 

CALL SINTRP(RS,CQ,FRENO,NCQ,lS, II I. I ORDER ,COD»NOPRNT ) 
DMSAVE=DM ( K ) 

GO TO AO 
601 CONTINUE 

ROAVG=(RO+RHO(K) )*.5 
VTR I L=A8S ( VV) 

RE=(R0AVG*VTRIL*2.*RA0(K) )/( AMUU*GC) 

C FRICTION COEFFICIENT FOR SMCCTH PIPES 

C ALLOW CHANGE OF RE RANGE UP TO 10 STEPS 

IF (COUNT .LE. 10) GO TO 70A 
C FORCE CALC TO REMAIN IN ONE RE RANGE 

I F ( ICONS -2)702,700,701 

7 OA I F ( RF .LT. 100000. .ANO. RE .GT. 1185.) GO TO 700 
IFIRE .LE. 1185.) GO TO 701 
C NIKURADSE EQN 

702 FF=(.0082+. 05 525/RE **(.237) ) 

ICCNS=1 
GO TO 705 

C EQN FOR BLASIUS RANGE 

700 FF=. 0791/RE**. 25 
ICCNS=2 

GO TO 705 

701 FF=16./RE 
I CCNS=3 

705 CONTINUE 

LOVRD=AL(K)/RAD(K) 

C CALCULATE (1/RO-l/RHO) FOR COMPRESSIBLE FLOW TUBE FLOW EQN 

ROCOMP=ABS( ( 1 ./RHO(K ) )-( 1 ,/RO ) ) 

C COMPRESSIBLE FRICTION FLOW EQN INVERTED TO SOLVE FOR A 

C VELOCITY DENSITY-COUPLE USED TO CALC MASS FLOW RATE 
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VV=SQR T( (ABS(P(K)-PM) / (FF*RCAVG*L0VRD+R0AVG**2*R0C0MP ) )*GC > 
C0UNT=C0UNT+1. 

C COO REPRESENTS THE FRICTION FACTOR FF WHEN TORE FLOW IS USED 

CQC=FF 

l F ( VV .NE. 0.) GO TO 1050 
VV=1. 

OM I K ) =0. 

GO TO 1060 

1050 IFIABSt VV/VTRIL-1. ) -.01 >600,600,601 
600 IFIVV-SS)630,630,6A0 

640 VV=SS 

C COO FOR M=1 IS ONLY APPROX I PA TEC BY .611 IT SHOULD BE 

C INVESTIGATED IF ANY AMOUNT OF CHOKED ORIFICE FLOW IS ANTICIPATED 

IFINOQ .EQ. 2 .OR. NOQ .EQ. 3> VV=AB*CQD 
630 IF(P(K)-PM>602,603,603 

602 VV=-VV 

603 0M(K)=R0AVG*PI*RA01K)**2*VV 
1060 DMTOT=DMTOT+OMlK) 

SUM=SUM+DM(K)*(CCP* TE) 

COW ( K ) =COD 
CPWIK)=CP 
REW ( K ) =RE 
PH 1 1 I K ) = PH I 
V VCUT ( K ) =VV 

C CALCULATE APPROXIMATION TO HF AT TRANSFER FOR PIPE FLOW 

C AL AT=L ATERAL AREA OF TUBE , QA 1 R=BTU/ SEC , TCONV=AVG TEMP OF AIR IN 

C TUBE, TUBE=AVG TUBE TEMPAVOL=VOL CF AIR IN EACH TUBE,TVOL= 

C VOL OF TUBE MATERIAL 

IFINOQ .FQ. 0) GO TO 305 

IFITCON .LE. 0. .OR. NOO .GE. 2> GO TO 305 
IFIRE .GT. 1185.) GO TO 310 

C CALC LAMINAR H, THRCON=THERM AL CONDUCTIVITY OF AIR 

HA = 2 • 1 82 *THRCON/R AD I K > 

GO TO 300 

C CALC TURBULENT H 

310 HA=ABSIVV)*R0AVG*CCP*C0D*.5 

300 CCNTINUE 

OTAIR=lQAIRIKJ*OELT)/<CCP*ROAVG*AVOLlK> > 

DTUBE=-IOAIRIK) *OELT) / 1 C P TUBE *RCTUBE*TV CL 1 K > > 

TME AN= l TE+TCQNV IK ) ) *.5 
TCCNVl K ) =TMF AN+OT AIR 
TUBEIK)=TUBE(K)+DTUBE 
C QAIR IS NEGATIVE CUT OF THE AIR 

QAIRIK)=HA*ALATIK) *1 TUBE IKJ-TMEAN) 

C CONVECTIUN HEAT TRANSFER FROM-TO AIR- TUBE 

IFtPIK) .GT. PM) QTOT = QTOT-*-OAIRIK) 

305 CCNTINUE 
500 CONTINUE 

C INTEGRATE EONS HERE ,SENO VALUES FOR CALC DERIVATIVES 

C INTO FUNCTION SUBROUTINE 

C ROD=OMTOT/VR 

C TMC =1 SUM-CV*TM*DMTOT*QTCT )/!CV*RO *VR) 

C 1 1) =DMTOT 
C I 2 ) = VR 
C I 3 ) =SUM 
CIA) =CV 
C I 5 ) =QW 
C I 5 ) = QTOT 
CALL CONSTSIC) 

20 CCNTINUE 
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J=GNE 

CALL INTEGT! J,RO,ROD,T,I SAVE.DELT) 

J=TWO 

CALL INTEGT! J,TM,TMD,T, I SAVE.DELT) 

I F ( I SAVE .EQ. 1) GO TO 20 
22222 CONTINUE 

C COMPUTE TOTAL MASS CHANGE IN MANIFOLD UP TO THIS STEP 

TOTMAS=TOTMAS+DMTOT*DELT 
I F ( TISO .LT. 0.) TM=ABS ( T I SO ) 

PM=RO*R*TM 

IF ( ( T-TPRNT ) .LT. PRNTFQ .AND. NCURVS .GT. 0 .AND. PLOTFQ .NE. 

1 PRNTFQ .AND. T .LT. TENOJ GO TO 10201 
IF ( ( T-TPRNT ) .LT. PRNTFQ .AND. T .LT. TEND) GO TO 21 
1 = 1+6 

IF! I .LT. 60) GO TO 505 
WRITEI6, 10100) 

WRITEI6, 10101) 

HRITEI6, 10106) 

WRITE(6, 10105) 

10100 FORMAT ( • 1 ALT I TUDE • , 2X , • VELOC I T Y • , 7X , • T I ME • , 3X, • MAN I FOLD PRESS* , 
IAX, ‘PRESSURE (PSF) AT EACH PORT* ,2AX, • ATMOS PRESS’ , AX,* QTOT* ) 

10101 F0RMAT('0*,2X,'00* ,8X,'MACH NO* , 8X ,* ALPHA* ,2X ,* MANI FOLD TEMP*,5X,* 
1PRESSURE COEFFICIENT AT EACH PCR T* , 1 8X , • ROLL RATE*) 

10106 FORMAT! *0 MASS FLOW RATE IN EACH TUBE • ,23X , 'REYNOLDS NO IN EACH T 
1UBE',27X, ’ INTERNAL SS* ,AX , • E X TERNAL SS* ) 

10105 FORMA T( ' 0 TEMPERATURE AT EACH PORT* ,26X,' PHASE ANGLE OF EACH PORT 
1* ,27X,* TOTAL MASS CHG IN MANIFCLD*) 

10103 FORMAT! *0 VELOCITY IN EACH TUBE • ,29X , • D ISCH ARGE COEFFICIENT AT EA 
1CH PORT' ) 

10107 FORMAT! *0 VELOCITY IN EACH TUBE* ,29X, • FRICTION FACTOR IN EACH TUB 
IE * »23X, • ODOT 1*,8X,*CD0T 2*) 

10108 FORMAT! *0 TCONV, MEAN GAS TEMPERATURE IN EACH TUBE* ,10X, • MEAN TUB 
IE TEMPERATURE OF EACH TUBE ' , 17X , • QCC T 3*,8X,*CD0T A* ) 

10109 FORMA T( ' 0 A VOL ' ,G 13 . 6 , • ALA T» ,G 1 3 . 6 , • TVCL* ,G1 3 . 6 ) 

IF! NOQ .EQ. 0 .OR. NOQ .EQ. 1) WRI TE ! 6 , 10 107 ) 

IF! NOQ .EO. 2) WR I TE ! 6 , 10103 ) 

1 = 10 . 

I F ! TCON .NE. 0. .AND. NOC .EQ. 1) WRI TE !6, 10108 ) 

IFtTCON .NE. 0. .AND. NOQ .EQ. 1) 1=1+1 

505 CONTINUE 

WRITE! 6, 1011 0)HO, V 1 , I , PM ,P ! 1 ) , P ! 2 ) , P! 3 ) , P ( A ) .PO.QTCT 
WRITE!6, 1010A)Q0, AA,ALPHA,TM , ( CPW < l Z ) , I Z= 1 » A ),OMEG 
WRITF(6,1010A)!DM(IZ)»IZ=1,A ) , (REW! IZ),IZ=1»A ),SS,AB 
WRITE16, 1010A) (TA( IZ) , I Z = 1,A ) , ! PHI III Z ) , l Z= 1 , A ),TOTMAS 
WRITEI6.1010A) ( VVOUT! IZ) , I Z = 1 , A ) , ! C C W ( I Z) , I Z = 1 , A ) 

I F I TCON .NE. 0. .AND. NOQ .EO. 1) WR ITE ! 6, 10 102 ) QA I R I 1 ) ,QAIRI 2 ) 

I F (TCON .NE. 0. .AND. NOQ .EQ. 1) WRI TE ( 6, 1 010A ) ( TCONV! IZ ) , I Z=1 , A 

l) , (TUBE! IZ) * I Z= 1 * A ) , CAIRO) ,QAIR!A) 

I F ( TCON .NE. 0. .AND. NOQ .EQ. 1) 1=1*1 

10102 FORMAT! * + * , 10AX , 2G1 3 . 6) 

1010A FORMAT! 1X.10G13.6) 

10110 FORMAT(IHO,10G13.6) 

I F ( NCURVS .EQ. 0) GO TO 90 

C VARIABLE PLOT POINT CAPABILITY IS CHECKED HERE 

10201 I F ( PLOTFQ .EO. PRNTFQ) GO TO 10202 
PLCTF=PL0TF+1. 

I F ( PLOTF .LT. PLOTFQ) GO TO 10203 
PLCTF=C. 

10202 CONTINUE 
PLCT ( 1 ) =T 
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PLOT ( 2 ) = PM 
PLC T ( 3 ) =P ( 1 ) 

PLC T ( 4) =P I 2 ) 

PLOT ( 5 ) = P I 3 ) 

PLOT ( 6) =P ( 4 ) 

PLOT ( 7 ) =PO 
PLOT ( 8 ) =HO 
PLOT ( 9 ) =V 1 
PLC T ( LO) =OTOT 
PLOT (11) =Q0 
PLOT ( 12 ) = AA 
PLOT ( 13)=ALPHA 
PLOT ( 14 ) = TM 
PLOT ( 15 ) = REW I 1 ) 

PLOT ( 16 ) = REW ( 2 ) 

PLOT ( 17 ) = REW ( 3) 

PLOT ( 18) =REW(4) 

PLCT(19)=CQW(1) 

PLOT ( 20 ) =CQW I 2 ) 

PLOT ( 21 ) =CQW ( 3) 

PLCT(22)=CQW(4) 

00 10200 I PL = l.NCURVS 
10200 APLOT ( IPL)=PLOT ( I PLOT ( I PL) ) 

CALL RJPLOT (APLOT *0) 

10203 IF ( ( T-TPRNT) -LT. PRNTFQ) GO TC 21 
90 CONTINUE 

I F ( T .LT. TEND) GO TO 35 

1 = 58 

I F ( PL TCON .EO. 1.) GO TO 10209 
CALL RJPLOT ( APLOT ,-l) 

10209 RETURN 
ENC 


/* 

//SYSLMOD ^ 00 UNIT=0ISK,V0L=SER=N2SCR6,0SN=N0JFLMFD.DISP= (NEW, KEEP) 
//LINK. OBJECT OD * 

.. DATA QMS iNi JLMFO 
.. DATA CMS.N.RJPLT 
.. DATA DMS * N» USUBS 
ENTRY MAIN 
NAME NOJFLMFDIR) 

/* 

.. DATA OMS . N « GC A19 
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SUBROUTINE CSUBP 

COMMON/ Z ILCHl/YXI 15 > , YPR IMEI 15» .ERCRIT ( 2 ) ,P< 41 , RADI4) , AL (4) . 
1ALT(50),VEL(50) ,TF(50),CQ< 10),FREN0( 10),CMACHl 10), ALFA! 10,15), 
2ALPHAAI 10) ,TALPH( 10) , PRESSF ( 6 ) ,TCCN, TTUBE, THICK, CPTUBE.ROTUBE, 

3F l (10, 15, 20) ,CPP( 10,15,20) ,T ,DELT , DELMAX, PM, VR, THETA, OMEGA, 

4THRC0N* PLTCON , EAZANG, ALPHA, CP, AA.TM, PRNTFQ, TEND, PHI , RT, INDPVR, 

5 JMAX , N0NL1 N *IERCRT»IXMAX*IYMAX»1ZNAX*KT»N0G» NCPTS, N CO, N ALPHA, 
6NPPTS*NCURV$» MODE ,IPSKIP,IPL0T(7) 

PRESSURE COEFFICIENT SUBROUTINE 

CP IS DETERMINED FROM M NO, ALPHA, AND PHI FROM EXPERIMENTAL CURVES 
THREE WAY LINEAR INTERPOLATI CN USED TO FIND CP 
PHI IS THE ROTATION PLUS PORT POSITION AND IS SET .LE. 180 DEG 
SINCE CP CURVES ARE CONSIOEREO SYMMETRIC 

AT IS EITHER MACH NO OR TIME DEPENDING ON THE VALUE OF 
PRESSF ( 6 ) , IE, PRESSF (6 )=0, 1 FOR MACH NO, =-1,2 FOR TIME 
AT = AA 

IF ( PRESSF ( 6 ) .EQ. -1. .OR. PRESSF I 6) .EQ. 2.) AT=T 
PI = 3. 14 159 
ALPHA I=ALPHA 
I XN I N= l 
I YM I N= 1 
l ZM I N= 1 
IX=1 
I Y = 1 
IZ=l 

M I N TRP = 0 

IF TABLE END POINTS ARE HIT MINTRP =2 STOPS DOUBLE INTERPOLATION 
I F ( AT .LT.CMACHI IXMAX) .AND. AT . GT . CM ACH I I XM IN ) ) GO TO 501 
I F ( AT .LE. CMACHI IXMIN) ) GO TO 504 
I X= I XMA X 

CMACHI IX) =CM ACH (IXMAX) 

M I NTRP=2 
GO TO 503 
04 IX= IXMIN 

CM ACH ( IX)=CMACHI IXMIN) 

MINTRP=2 - 
GO TO 503 

01 CONTINUE 

IF (AT .GE. CMACH(IX) .AND. AT .LT. CMACHIIX+D) GO TO 502 

IX=IX+1 

GO TO 501 

02 CONTINUE 

INTERPOLATE ON MACH NO, MINTRP=1 TRIGGERS SECOND INTERPOLATION 
405 IFIMINTRP .EQ. 1) IX=IX+1 
503 CONTINUE 

I F ( ALPHA 1 . L T. ALF A ( I X , I YM AX ) .ANO. ALPHAl .GT. AL FA ( IX , IYM IN ) ) 

1 GC TO 511 

IF( ALPHAl .LE. ALFA ( I X , I YMI N ) ) GO TO 514 
I Y= I YMAX 

ALFA(IX,IY)= ALFAIIX.IYMAXI 
GO TO 512 
514 IY=IYMIN 

ALFA ( I X , I Y) = ALFA ( I X , l YMI N) 

GO TO 512 

511 IF( ALPHAl. GE. ALFA(IX.IY) .AND. ALPHAl .LT. ALFA ( I X, I Y+ 1 ) » GO T0512 
I Y= I Y+ 1 

GO TO 511 

512 CONTINUE 

C PH I = VEHICLE ROLL ANGLE 

PHI=OMEGA*T*THETA 
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RESET=PHI/(2.*Pl ) 

IRESET=RESET 

PHI=(RESET-l RESET f *2.*PI 

C FOR NONSYMMETRIC CP DATA SKIP ASSUMPTION OF SYMMETRY 

I F ( PHI .GT. PI .AND. PRESSF ( 6 ) .LT. 1.) PHI = 2.*PI-PHI 
C PHI IN DEG = PHIIR AO 1*57.2961 DEG/ RAD ) 

PHI=PHI*57.296 

540 I F ( PH I .GE. FIIIX.IY.IZ) .ANO. PHI .IE. FIUX.IY.U+l)) GO TO 550 

IZ=IZ+1 
GO TO 540 
550 CONTINUE 

C TEMPORARY PROCEEDURE FOR SOLVING END POINTS 

I F ( I Y .LT. IYMAX) GO TO 555 
I Y= I YM AX- 1 

555 I F ( I Z .LT. IZMAX) GO TO 560 
I Z = I Z MAX- 1 
560 CONTINUE 

D I V 1=F I ( IX, IY + 1, IZ + 1 )-FI ( IX,IY+1,IZ) 

0IV2=FI I IX.IY ,IZ+1)-FI (IX, IY ,IZ) 

SLCPE 1= ( CPP (IX, IY+1 , IZ+1 l-CPPI IX, IY+1, IZ))/DIV1 
SL0PE2= ( CPP ( IX, IY, IZ+l )-CPPI IX, IY, I Z ) )/DIV2 
CP1= SLOPE 1*( PH l -FI { IX,1 Y,IZ) )+CPP( IX,IY+l,IZ) 

CP2= SL0PF2* (PH I-FI ( IX.IY, IZ) ) +CPP ( I X , IY , I Z ) 

DIFF=(CP1-CP2)/(ALFAI IX, IY+1 )-ALFA( IX.IY)) 

563 ALPH= ALPHA1 

I F ( ALPH .LF. ALFA! IX, IY+1) ) GO TO 570 
ALPH=ALFA( IX.IY+l ) 

570 CP= OIFF*( ALPH-ALFAIIX, IY) )+CP2 
IF(MINTRP .GT. 0) GO TO 400 
M I NTRP= 1 
CPI=CP 
GO TO 405 

400 IFIMINTRP .EO. 2) GO TO 410 
C CP VALUE 8FTWFEN MACH NUMBERS 

CP=( ICP-CPI )*( AT-CMACHI IX-1 ) ) UtCMACHI IX)-CMACH( I X— 1 ) > +CPI 
410 CONTINUE 

C FAZANG IN RAD =2.*PI/KT) WHERE KT = NO CF PORTS 

THETA=THETA+FAZANG 
RETURN 
END 
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SUBROUTINE RE AD IN < D I CT, B ,ND ICT ) 

REAL*8 DICT(NDICT),A(100), ARRAY 
INTEGER OPTION 

D l MENSI ON BINOICT), I NPUT ( 100 ) ,FMT (20 ) , FMT A (20 ) , 

1FMTWI20) ,FMTH0L(20) 

C INITIALIZE REAO AND WRITE FORMATS 

DATA FMT/' ( 10G' , • 8.5) • ,18*' •/ 

DATA FMTA/' ( 10G* » ' 8.5) • ,18*' •/ 

OATA FMTW/' ( 1CG' 13.6* ) ',17*' •/ 

C NOICT IS THE NUMBER OF DICTIONARY ENTRIES, IE HIGHEST INDEX USED 

DO 1000 1=1,100 

A ( I ) =0. 

INPUT! I ) =0 
1000 CONTINUE 

1 FORMAT ( 1615) 

2 FORMAT ( 10A8 ) 

3 FORMAT ( 20A4 ) 

FORMAT! IHO, 10A8) 

FORMAT !A8, I 7, 1315, !/,16I5) ) 

0 FORMAT ( • 0 THE NAME ',A8,'IS NOT IN THE DICTIONARY, IT IS MISSPELL 
1ED OR MISPLACED IN THE FIELD. PLEASE CORRECT AND RESUBMIT PROGRAM') 

ISTOP IS USED TO STOP THIS PROGRAM WHEN INPUT NAMES ARE 
UNRECOGNIZABLE , ALL NAMES ARE PROCESSED BEFORE STOPPING 
I STOP=0 

READ IN THE NUMBER OF VARIABLE NAMES OR NO OF VARIABLES IN AN ARRAY 

01 CONTINUE 
NEWFMT=0 

READ! 5, 1 )OPTION,NOVAR,NEWFMT 

NEWFMT SETS UP OUTPUT FORMATS 

0 - NO FMT , 1- WRITE FM T, 2-HOLLER I TH FMT ,3-BOTH 1C2 
I F ( NEWFMT .EO. 2 .OR. NEWFMT .E0.53) READI5.3) FMTHOL 
I F ! NEWFMT .EO. 2 .CR. NEWFMT .EQ.53) WRITE!6, FMTHOL ) 

I F ! NEWFMT .EQ.51 .OR. NEwFMT .EQ.53) READ!5,3) FMTW 

OPTIONS- 0= READ VAR NAME ( S ) -RETURN, 1=0+READ ANOTHER OPTION 
AFTER THE DATA, 2=READ IN AN ARRAY NA ME-RE TURN , 3=2+R E AO ANOTHER 
OPTION, 6=READ IN AN ARRAY NAME ♦ SPECIFIC ARRAY LOCATIONS, 7= . 
6+ READ IN ANOTHER OPTION, A PREFIXES ANY PREVIOUS OPTION NO 
WHEN IT IS DESIRED TO WRITE THAT INPUT DATA, IE 40,46, ETC, 

5 PREFIXES ANY PREVIOUS OPTION TC INPUT A NEW DATA FORMAT 
IE, RE-DEFINE A PREVIOUSLY OFFINED DATA FORMAT, 50,56,540,543 
NOTE OPTION 4 MAY NOT PRECEDE OPTION 5 IE, 450 IS INVALID 
I SK I P=5 
MINUS=50 

I F ( OPT I ON .LT. 50) GO TO 77 
I F ! OPT I ON .GF. 500) MINUS = 500 
C SET OPTION 

OPT I ON=OPTI ON-MINUS 
I SK I P=0 

77 CONTINUE 

I F (OPT I ON .EQ. 2 .OR. OPTION .EQ. 3) GO TO 102 

IF! OPTION .EQ. 42 .OR. OPTION .EQ. 43) GO TO 102 

I F ( OPT I ON .EQ. 6 .OR. OPTION .EQ. 7 .OR. OPTION .EQ.46 .OR. 

1 OPTION .EQ. 47) GO TO 103 
C READ VARIABLE NAMES INTO THE A ARRAY BY A FORMAT 
REAO! 5,2) (A(I),I=1, NOVAR ) 

DO 5 1=1, NOVAR 
00 6 K= 1 , ND ICT 

IF! 0 IC T { K ) .EQ. AID) INPUT!I)=K 
IF! DICT(K) .EO. AID) GO TO 5 
6 CONTINUE 
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C IF THIS WRITE IS IMPLEMENTED AN UNRECGONI ZABLE NAME HAS 

C BEEN FOUND, NAME PROCESSING WILL CONTINUE UNTIL END OF NAMES 

WRITEI6.10) AID 
I STOP=l 
INPUTII )=1 
5 CCNTINUE 

C READ THE FORMAT FOR THE VARIABLES TO BE READIN 

IFIISKIP .NE. 5) READ (5,3) FMT 
READ( 5, FMT ) ( B( INPUT! I ) ) ,1=1, NO VAR) 

I F ( OPT I ON .EO. 1) GO TO 101 

I F ( OPT I ON .EQ. 40)WRITE(6,FMTW)( B UNPUT (I) ) , I = 1 , NOVAR > 

I F (OPT ION .EO. 41 )WR!TE(6,FMTWM B ( INPUT ( I )), 1= 1, NOVAR ) 

I F ( OPTION . EQ.41 ) GO TO 101 
IF( ISTOP .EQ. 1 ) GO TO 75 
RETURN 

C NOARRY IS THE NO OF ARRAY VALUES TO BE READ IN STARTING AT 
C ARRAY LOCATION ISTART 

102 NOARRY=NOVAR 

RE AD ( 5,9) ARRAY, ISTART 
I F ( I START .LT. 1) ISTART = 1 
DO 50 K= 1 , NO ICT 

I F ( ARRAY .EQ. DICT(K)) GO TO 55 
50 CCNTINUE 

C IF THIS WRITE IS IMPLEMENTED AN UNRECGONI ZABLE NAME HAS 

C BEEN FOUND, NAME PROCESSING WILL CONTINUE UNTIL END OF NAMES 

WR I TE ( 6 , 10 ) ARRAY 
I ST0P=1 
K=1 

55 KA = KUSTART-1 

NPLACE=KA+N0ARRY-1 

IFIISKIP .NE. 5) RE AD ( 5 , 3 ) FMT A 

READ( 5.FMTA) ( B (I) , I =K A ,NPLACE ) 

C LOAD B ARRAY WITH ARRAY VALUES FOR TRANSMISSION TO USERS PROG 

I F ( OPT I ON .EQ. 3) GO TO 101 

I F ( OPT I ON .EO. 42) WR I T E ( 6 , FMTW ) ( B (I) , I = KA, NPLACE ) 

I F ( OP TI ON .EO. 43) WR I TE ( 6 ,F MTW ) ( B (I) , I = K A, NPLACE ) 

I F ( OPT I ON .EO. 43) GO TU 101 
I F ( I STOP .EQ. 1) GO TO 75 
RETURN 

103 CCNTINUE 

RE AO ( 5 , 9) ARRAY, ( INPUTU) , 1 = 1, NOVAR) 

DO 60 1=1 , ND I CT 

I F ( ARRAY .EQ. DICTUM GO TO 61 

60 CCNTINUE 

C IF THIS WRITE IS IMPLEMENTED AN UNRECGON I Z ABLE NAME HAS 

C BEEN FOUND, NAME PROCESSING WILL CONTINUE UNTIL END OF NAMES 

WRITE(6,10) ARRAY 
I STQP=1 
1 = 1 

61 CCNTINUE 

I SA VE=I- l 

DO 62 1=1, NOVAR 

INPUT! I ) = INPUT I I ) ♦ ISAVE 

62 CCNTINUE 

IFIISKIP .NE. 5) RE AD ( 5 ,3) FMTA 

READ ( 5, FMT A ) ( B( INPUT ( I ) ) , 1=1, NOVAR) 

C LOAD B ARRAY 

I F ( OPT I ON .EO. 46 .OR. OPTION .EQ. 47) WR ITE ( 6, FMTW ) ( B I I NPUT ( I ) ) , 
1 1=1, NOVAR) 

I F ( OPT I ON .EO. 7 .OR. OPTION .EO. 47) GO TO 101 

75 l F ( ISTOP .EQ. 1) STOP 

RETURN 
ENO 
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SUBROUTINE SINTRP (A,B,C , I M AX, I S , I , I ORDER, D .NOPRNT ) 

IMPLICIT REAL*8(A-H,C-Z) 

R£AL*A A, B,C,D 
DIMENSION B< ICO), Cl 100) 

DIMENSION XI 100) ,Y< 100) ,NUM< 20),CEN0M( 20, 20),U(20) 

REAL NUM.NUMT 
X I = A 

22223 CCNTINUE 

DO 60 K=1,IMAX 
Y I K ) = B I K ) 

X<K)=CIK) 

60 CCNTINUE 

I F ( I ORDER .LE. 1) I0RDER=2 
C SEARCH AND INTERPOLATE ROUTINE 

C IS=0 FOR BOTH SEARCH AND INTFRP 

C I S= 1 FOR SEARCH ONLY 

C I S=2 FOR INTERP ONLY 

C N0PRNT=3 SUPRESSES THE WRITTING OF THE 

C DIAGONISTIC MESSAGE , NOPRNT RETURNS A VALUE OF -1,0, OR 1 DEPENDING 
C ON WHETHER THE INTERPOLATION TCCK PLACE OUTSIDE THE TABLE .LOWER 
C (-1 ) ,UPPEREND< 1) ,0R INSIDE THE TABLE ( 0 ) 




IFIIS 

.EQ. 

2) 

GO TO 

20 



1 = 1 







IF! XI 

.LT. 

XI 

IMAX) 

. ANO. XI .GT. XI l) ) GO TO 10 



IF I XI 

• GT. 

XI 

1) ) GO 

TO IS 



1 = 1 







GO TO 

2C 




1 

5 

l = I MAX 







GO TO 

20 




1 

0 

I F I X I 

.LE. 

XI 

1 + 1) . 

AND. XI .GT.XI I ) ) GO TO 20 



1*1*1 







GO TO 

10 




2 

0 

CCNTINUE 






IFIIS 

.EQ. 

1) 

RETURN 


C TABLE WITHIN A TABLE INDEX LOCATION 

I T= 1 

IFIIS .EQ. 2 .AND. I .GT.1)IT=I 
IFIIORDER .GT. IMAX) IOROER=IMAX 
IORDRS=IORDER 

22 I F ( I I + 1 ORDER- 1 ) .LE. IMAX) GO TO 21 
IGRDER=IQRDER-l 
IFIIORDER .NE. 1» GO TO 22 
A6 IF (NOPRNT .EQ. 3) GO TO AS 

WRITE I 6, 101 ) Y( IT) ,X< IT) ,Y( IMAX) , XII MAX) ,Xl, IMAX 
101 FORMAT I IHO, • AN ERRONEOUS VALUE HAS BEEN SUPPLIED TO SINTRP, IE XI 

1 IS OUTSIDE THE TABLE BEING USED , Yl = ' ,015 . 8 , * Xl=’ ,C15.8 ,/ , * YHAX= 

2 •»Dl5.8,*XMAX=',Di5.8, , XI=* ,D15.8,*IMAX=*»I5»/»* THE APPROPRIATE 
3ENC POINT IS USED FOR THE INTERPOLATED VALUE*) 

A 5 CCNTINUE 

I I = I MAX 

I F I X I .LT. X I IMAX ) ) 1 1 = 1 
UX = Y I I I ) 

IQROER= IORDRS 
NOPRNT *— 1 

IF! II .EQ. IMAX) N0PRNT=1 

D=UX 

RETURN 

21 I F I X I .LT. X I IT ) ) GO TO A6 

C SAVE THE ORIGINAL INDEX I TO RETURN TO CALLING PROGRAM 

I RETRN= I 
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THE INDEX I REPRESENTS THAT POINT IN THE TABLE WHERE THE 
POLYNOMIAL WAS STARTED 
I F ( IORDER .NE. IORDRSI I = I - I IORDRS- I ORDER ) 

NTH ORDER INTERPOLATOR 
0 CONTINUE 

RETURN IORDER TO ORIGINAL VALUE 
I OROER= I ORDRS 
NOPRNT=C 
I SAVE =1 

DC 25 11=1, IORDER 
I = I S AVE+ I I-l 
NUM ( l I ) = XI-X( I) 

DC 25 LL=1, IORDER 
L = LL + 1 SAVE- 1 
DENUHIII ,LL)=X( I ) -X ( L ) 

C SET DIAGONAL =1 SO THAT DIVISION BY 0 DOES NOT OCCUR 
I F ( I .EQ. L) DENOM ( 1 1 , LL ) = 1 . 

25 CONTINUE 

DO 26 1 = 1, IORDER 
TERM =1. 

DO 27 L=l, IORDER 
I I =L 

I F ( I .NE. L) GG TO 30 
NUMT = NUMI 1 1 ) 

NUM ( I I ) = I . 

30 TERM = TEKM*(NUM( II )/DENOM(I,LU 

I F ( I .EQ. L) NUM ( II) =NUH T 
27 CCNTINUE 

U « I ) =TERM 

26 CCNTINUE 
1 = 1 SAVE 

C VALUE OF I NOE P VAR AT XI 
UX = 0. 

00 AO K= l, I ORDER 
UX=UX+U( K) *Y( I ) 

1 = 1+1 - 

AO CONTINUE 

1= I RETRN 
0 = L X 

22221 CCNTINUE 
RETURN 
END 
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SUBROUTINE AT62 <ZFT,ANS) 

REAL MOLWT.LOGPR 

DIMENSION H ( 23 ) » TBASE I 22 ) , TGRAD < 22 > ,PBASEI 22 ) .MCLWTI23) ,ANS<4) 
DATA H/O., 11., 20., 32., 47. , 52., 6 1., 79., 88. 743, 90., 100., 1 10., 120., 1 
1 50. » 160 . * 1 TO. » 190 . * 230. , 300. *400 .*500. ,600., 700. /,T3ASE/288. 15*216 
2.65,216.65,228.65,270.65,270.65,252.65, 180.65, 180. 65, 1 80. 65, 210. 65 
3,260.65,360.65,960.65,1110.65,1210.65,1350.65, 1550.65,1830.65,2160 
4.65,2620.65,2590.65/ ,TGRAD/— 6.5,0. ,1», 2.8,0.,— 2.,— 4«,0.,0.,3»,5«,1 
50., 20. ,15., 10. ,7. ,5. ,4. ,3. 3, 2. 6, 1.7,1. 1/»PBASE/ 1 . , 2. 2336 IE— 01, 5.40 
6328E-02,8.566636-03,1.09455E-C3,5. 8 2 28 9 E -04, 1.79718E-04, 1.0241 E-05 
7.1.6223E-06, 1 .622 3E-06 ,2 .9681E-07 , 7. 2582E-08 , 2.4887E-08, 4.9955E-09 
8,3.6460E-09,2.7561E— 09, 1 . 6632E-09, 6. 8694E- 1C , 1 . 8592E-10, 3.9777E- 11 
9, 1. C814E- 11 , 3. 405E- 12/ ,M0LWT/1C*28. 9644 ,28. 88,28.56 ,28.07,26 .92 , 26 
A. 66, 26. 40, 25. 85, 24. 70, 22. 66, 19.94,17.94, 16.84, 16. 17/ , RE/6355.63/ ,C 
B/34.1628/,CRH0/.3236458E-03/,PSFA/2116.22/ 

IF (ZFT.LT.O.) ZFT=0. 

IFIZFT .GT. 2.E06) ZFT=2.E06 
Z=.3048E-03*ZFT 
IF (Z.GE.90. ) GO TO 5 
Z=Z/ ( 1. + ( Z/RE) I 
DO 1 J=1 ,9 

IF (Z.GE.HI J). AND.Z.LE.H! J*l) ) GO TO 2 

1 CONTINUE 

2 CONTINUE 

TKELV= TBASE! J) * TGRAD ( J) MZ-H! J) ) 

IF (ABSITGRAD(J) ) .L T . .5 ) GO TO 3 

PSF=PBASE!J)*PSEA*( ( TBASE ( J ) /TKELV ) ** (C/TGR AD ( J> ) ) 

GO TO 4 

3 CONTINUE 

PSF=PBASF< JJ*PSEA*EXP!-C*!Z-H< J> I /TBASE! J) ) 

4 CONTINUE 
$LGFT3=CRHU#PSF/TKELV 
VS0UND=S0RT< 4325.746*TKELV) 

GO TO 8 

5 CONTINUE 
00 6 K= 10, 22 

IF (Z.GE.H(K). AND.Z.LE.H(K-UJ) GO TO 7 
CONTINUE 
CONTINUE 

TM=TBASE!KM-TGRAD!K) *< Z— H< K) ) 

A* 1 .♦ ( H ( K ) /RE ) 

B=TBASE!K)/!TGRAD(K)*RE) 

X=!Z-H(K) I/RE 

LOGPR=(-C/TGRAD(K))>M l./C A-B) )*U1./(A+X ) ) - { 1 ./ A ) + ( 1 ./ ( A-B ) )* ALOG 

1 ( ( A* ^ 8+X ) ) / ( B*( A+X) ) ) ) 

PSF=PBASE(K)*PSEA*EXP(LOGPRI 
SLGF T3=CRH0*PSF / TM 

DMOLWT= ( (MOLWTIK+1) -MOLWT < K ) J / ( H !K + 1 » -H ! K ) ) )*(Z-H(K) > 

TKELV= ( ( MOLWT! K) +DMCLWT 1/28.9644 )*TM 
VS0UND=894.50 
8 CONTINUE 

ANS l 1 )=SLGFT3 
ANS ( 2 ) =PSF 
ANS! 3I=TKELV 
ANS! 4) =VS0UND 
RETURN 
END 
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SUBROUTINE START 
IMPLICIT REAL*8 (A-H.O-Z) 

RE AL*A WW,FDT,TT,DELTT,YX(15),XX,DELMX,ERCRIT<2) ,YDOT ! 15), Cl 151 
COMMON/ZILCH/ FOOT 1 10 , 15 ) , Wt 10 . 15 ) ,CST ( 15 I , YY I 1 0, 15 ) . I ( 15 ) ,T ( 10 ) , 

1 ZIP! 15) .FOOTS! 10, 15), PCER ,DE LT ,0E LMAX ,DEL , S AVE , S I , SAVE1 . SII, 

2 I SKIP! 15) , K , N, J, JMAX , I STEP, I DUBLE , ICHECK, NOOBL E , NONL IN, ITERAT 

C w = DEP VAR,T= INOEP VAR, DELT=CHANGE IN INOEP VAR 

C JMAX= NO OF EONS ,K,N,I, ARE FREE INOICIES 

C DELMAX = MAXIMUM ALLOWABLE TIME STEP-USER INPUT 

C YX= INITIAL VALUES OF THE INDEPENDENT VARIABLE 
C YPR l M= INITIAL VALUES OF YOOT FOR THE DEPENDENT VARIABLE 
C OELTT DELTA T .STEP SIZE INITIAL VALUE 
C YD0T1 I ) = I N I T I AL VALUES OF THE DERIVATIVE 

C ERCRIT(l)= FRROR CRITERIA FOR ST A8 I L l T Y , E RCR I T ( 2 ) = PRED-CORR CRITERIA 
C NONL l N= INDICATES WHETHER OR NOT ITERATION TECHNIQUE IS TO BE USED 

C I ERCRT== INDICATES WHETHER OR NOT THE DEFAULT OPTION FOR THE ERROR 

C CRITERIA IS TO BE USED, I ERCRT=0 ,1,2,3 - DEFAULT OPTION 

C , SII CHANGED, PCER CHANGED, BOTH SII AND PCER CHANGED 

ENTRY STAR Tl ( YX, YDO T , XX ,0E L T T ,DE L MX , E RC R I T , J AX ,NCNLN , I ERCRT ) 
JMAX= JAX 
. NCNL I N=NONLN 
DELT=DELTT 
OELMAX=DELMX 
Til) = XX 

C DEFAULT OPTION FOR ERROR CRITERIA 
SII=.01 


5 


20 


PCER=.Cl 
IF! I ERCRT .EC. 
IF! I E RCRT .EQ. 
I F ( I ERCRT .EQ. 
IF! I ERCR T .NE. 
S I I = ERCR I T ( 1 ) 
PCER=ERCR I T ( 2 ) 
CONTINUE 
DO 20 1 = 1, JMAX 
Wt 1, I ) =YX ( I ) 

I F ( NONL I N .GT. 
ISKIP! I )=1 
ZIP! I ) =0. 
CONTINUE 
DEL=OELT 
SAVE1=0. 


0) GO TO 5 

1) SII =FRCR I Til) 

2 ) PCER=ERCRIT<2) 

3) GO TO 5 


0) FDOT ( 1 , I ) =YDOT ( I) 


K= l 

I DUBI, E= 1 
NCCBLE=0 


I S TEP= 1 
I SAVF=SAVE 1 
ICHECK=0 
RE L T I M= T ( 1) 

RETURN 

ENTRY CONSTS(C) 

C CST TRANSFERS CONSTS UR VAR TO THE DIFF EQN WHILE 
DC 30 I C= 1 , 1 5 
3 3 CST ( I C » =C ( I C ) 

RETURN 

ENTRY INTEGT! JJ.WW.FDT, TT , ISAVE, OELTT ) 


INTEG IS IN PROCESS 


J = JJ 

10 CALL INTGRT 
I S AVE = SAVE 1 


NN = N 

IF! J .LT. JMAX .AND. ISKIP(JMAX) .GT . A) NN=N+1 
WW=W(NN,J) 

FD T = FDOT I NN , J ) 

TT=T 1 NN ) 

DELTT=DELT 
22221 CONTINUE 
RETURN 
END 
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SUBROUTINE FUNCT 
IMPLICIT RE AL*8 (A-H,0-Z) 

COMMON/ZILCH/ FOOT 1 10, 1 5 > , M( 10, 1 5 » ,CST ( 15 ) , YY( 10, 15) ,Z 1 15 > . T ( 10 ) . 

1 ZIP! 15) , FOOTS! 10,15) ,PCER ,DELT ,DELMAX .DEL ,S AVE , SI ,S AVE1 , SII, 

2 ISKIPJ 15) ,K,N,J,JMAX,ISTEP, IDUBLE , ICHECK,NODBLE,NONL IN, ITERAT 
DIMENSION Y( 20) ,YOOT ( 20 ) 

J^EQNS ARE 0F r TH E FORM YOOT=F(X,Y), WHERE X IS THE INDEP VAR AND 
X=T ( K ) 

DO 50 1=1 , JMAX 
Y( I) = W(K, I ) 

I F ( K .LT. 8) Y( l ) =W( 1, I ) 

YOCTU ) = FOOT!K, I ) 

IF(K .LT. 8) YOOTI I )=FOOT( 1, I ) 

CONTINUE 

GO TO tl,2,3,4,5,6,7,8,9,10),J 
CONTINUE 

PLACE THE 1ST EQN ,YDOT=F!X,Y) IMMEDIATELY FOLLOWIN THIS CARD 
YDOTI 1)=CST(1)/CST(2) 

FDCT!K, J)=YDOT! J) 

RETURN 

CONTINUE 

PLACE THE 2ND EON ,YOCT=F(X,Y) IMMEDIATELY FOLLOWIN THIS CARD 
YD0T!2)=<CSTI3)-CST!4)*Y!2)*CST(i)+ CST { 5 ) )/ !CST!4)*Y! 1)*CST!2)) 
FDOT I K, J )=YDOT( J) 

RETURN 

CONTINUE 

PLACE THE 3RD EON ,YDOT=FIX,Y) IMMEDIATELY FOLLOWIN THIS CARD 
FOOT ( K, J ) =YDOT ( J ) 

RETURN 

CONTINUE 

PLACE THE 4TH EQN ,Y00T=F(X,Y) IMMEDIATELY FOLLOWING THIS CARD 
FDOT ( K, J )= YDOTI J ) 

RETURN 

CONTINUE 

PLACE THE 5TH EON ,YDOT = F(X,Y.) IMMEDIATELY FOLLOWING THIS CARD 
FOOT ( K, J ) =YDOT I J ) 

RETURN 

CONTINUE 

PLACE THE 6TH EQN ,YDOT=F(X,Y) IMMEDIATELY FOLLOWING THIS CARO 
FDOT I K , J ) =YOOT I J ) 

RETURN 

CONTINUE 

PLACE THE 7TH EQN ,YDOT=F(X,Y) IMMEDIATELY FOLLOWING THIS CARO 
FOOT ( K, J ) =YOOT I J ) 

RETURN 

CONTINUE 

PLACE THE 8TH EQN , YDCT=F I X , Y ) IMMEDIATELY FOLLOWING THIS CARD 
FDOT I K, J )=YDOT( J) 

RETURN 

CONTINUE 

PLACE THE 9TH EQN ,YDOT=FIX,Y) IMMEDIATELY FOLLOWING THIS CARD 
FOOT I K , J ) =YDOT I J ) 

RETURN 

CONTINUE 

PLACE THE 1 OTH EQN ,YDOT=F(X,Y) IMMEDIATELY FOLLOWING THIS CARD 
FDOT I K, J ) =YDOT I J ) 

RETURN 

END 
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SUBROUTINE ERROR 
IMPLICIT REAL*8 IA-H.O-Z) 

COMMON/ZILCH/ FOOT 1 10,15 ),W( 10,15 ),CST< 15), YYI 10*15). Z( 15 It Tl 10), 

1 ZIP! 15) , FOOTS! 10,15) ,PCER ,0ELT ,DELM AX, DEL .SAVE, S I , SAVE 1 , SII, 

2 ISKIP! 15) ,K,N,J«JMAX,ISTEP, IDUBLE, I CHECK, N0DBLE , N0NL IN, ITERAT 
0 I MENSI ON B! 1C) , AKI 10) 

DIMENSION Y ( 20 ) , YDOT ( 20 ) 

C DFDY REPRESENTS THE PARTIAL OF YOOT WITH RESPECT TO Y 

C THE EQNS ARE OF THE FORM DF0Y=F(X,Y), WHERE X IS THE INDEP VAR AND 

C Y IS THE DEP VAR 

C YDOT REPRESENTS THE DERIVATIVE CF THE FUNCTION BEING WORKED ON 
X=T IN+1 ) 

DO 50 1=1, JMAX 
Y(I)=W(N+l»t) 

50 YDCT ( I ) =F0CT ( N+ 1,1) 

C AK IS A CHECK FOR STABILITY, B CHECKS FOR TRUNCATION ERROR 

C AK=ABS(D(F) /DY) *OELT 

120 GO TO (1, 2, 3, 4, 5, 6, 7, 3, 9, 10), J 

1 CONTINUE 

C PLACE THE 1ST EON ,DFDY = F ( X, Y ) IMMEDIATELY FOLLOWIN THIS CARD 
DF0Y=0. 

AKIJ)=DABS(DFDY*DELT) 

GO TO 100 

2 CONTINUE 

C PLACE THE 2ND EON ,DF0Y=F(X,Y) IMMEDIATELY FOLLOWIN THIS CARD 
0FDY=0. 

AK! J)=DABS!DFDY*OELT) 

GO TO 100 
CONTINUE 

PLACE THE 3RD EQN ,DFOY=F!X,Y) IMMEDIATELY FGLLOW IN THIS CARD 

DFDY=!CSTl 1 ) +CST < 2 ) *DCOS < 4*Y < 4 ) ) ) **2*Y ! 3 ) - < C ST < 1 ) *CST ! 2 ) *DCOS ! 4. * 
1 Y ( 4 ) U*ICST!3)+CST(4)*DS!N(4.*Y!4) ) ) 

AK! J)=DABS(DFDY*OELT) 

GO TO 100 
CONTINUE 

PLACE THE 4TH EQN , DFDY=F ( X, Y ) IMMEDIATELY FOLLOWIN THIS CARD 
DFO Y=0. 

AK! J)=DABS!DFDY*OELT> 

GO TO 100 
CONTINUE 

PLACE THE 5TH EON ,DFOY=F(X,Y) IMMEDIATELY FOLLOWIN THIS CARD 
AK! J1=DABSIDFDY*DELTJ 
GO TO 100 
CONTINUE 

PLACE THE 6TH EQN ,DFDY=F (X,Y ) IMMEDIATELY FOLLOWIN THIS CARD 
AKI J)=DABS!DFDY*OELT) 

GO TO 100 
CONTINUE 

PLACE THE 7TH EON ,DFDY=F < X, Y ) IMMEDIATELY FOLLOWING THIS CARD 
AKI J)=DABSIDFOY*OELT) 

GO TO 100 
CONTINUE 

PLACE THE 8TH EQN , DFDY=F ! X, Y ) IMMEDIATELY FOLLOWING THIS CARD 
AK I J ) =DABS I DFDY+DELT ) 

GO TO 100 
CONTINUE 

PLACE THE 9TH EQN ,DFDY=F IX, Y ) IMMEDIATELY FOLLOWING THIS CARD 
AK ! J ) =DABS( DFDY *DELT ) 

GO TO 100 
10 CONTINUE 
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C PLACE THE 10TH fiON ,DFDY=F(X,Y1 IMMEDIATELY FOLLOWING THIS CARD 
AK(J)=DABS( DFDY#OELT ) 

100 CCNTINUE 

C Z< J) = PREOI J) /CQRRI J) 

B I J ) =DABS ( Z ( J ) ) 

56 IF ( J .LT. JMAXJ RETURN 

SAVE=B ( 1 ) 

. SI=AK( l) 

C DETERMINATION OF LARGEST ERROR TERM DETERMINES STEP SIZE 
JAX= JHAX-1 

I F ( JAX .EQ. 0) GO TO 205 

DO 35 L=l» JAX 

IF ( SAVE-B ( L + l ) ) 40.36, 36 

40 SAVE=B( L+l ) 

36 I F ( SI-AK1 L+l ) )4l,35,35 

41 SI=AK(L*l) 

35 CCNTINUE 

205 CCNTINUE 

SAVE= (DABS! SAVE I — 1» ) 

RETURN 

ENO 
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SUBROUTINE DERI V5 I X ,ZPC ,Vi .FOOT ,T , J, N) 

IMPLICIT REAL*8 (A-H.O-Z) 

EXPIQQQI-DEXPIQQO) 

ALOGI QQQ)=DLOG( QOQ ) 

ABS(QQQ)=DABS(QQO) 

DERIV ( A,B,C.D,E,F,G,H,0,P)=A/F*B/G*C/H-»D/0*E/P 
DIMENSION W( 10,15) .FOOT < 10,15) ,Y( 15),T( 10) 

DIMENSION DIVI20) 

DO 25 1=1,10 
25 DIVI I ) = 1. 

C IF THE DIVIDED DIFFERENCE EONS ARE NOT USED FOR THE 5TH DERIV 

C OF THE JTH EON SET DIV(J)=0. 

IF ( D I V ( J ) .EO. 0. ) GO TO 30 
FN*FDOT C N, J ) 

FN1=FD0T ( N-l , J) 

FN2=FD0T( N-2 , J) 

FN3=FD0T (N-3, J) 

FN4=FD0T (N-4, J) 

DENOM= (TIN)-TCN-l ) ) *( T I N )-T ( N-2 ) ) * I T IN )-T € N-3) ) *( T(N)-T(N-4) ) 
DEN0M1 = (TIN-U-T(N) )«( TIN-1) -TIN-2 ) )*(T (N-l )-T (N-3))* IT (N-l )- 
1 TIN-4)) 

0EN0M2=( TIN-2) -TIN) )♦ I T I N-2)-T (N-l ) )*( T IN-2)-T I N-3) )*(T( N-2)- 
1 TIN-4)) 

DEN0M3=( TIN-3 )-T< N) ) * IT IN-3 )-T (N-l ) ) *( T (N-3 )-T I N-2 ) )*(TIN-3)- 
1 TIN-4)) 

DEN0M4=(T(N-4)— TIN) ) * ( T I N-4 )-T (N-l ) ) * I T I N-4) -T I N-2 ) ) * I T I N-4 J - 
1 TIN-3)) 

IFIDENOM .NE. 0. . AND.DEN0M1 .NE. 0. .AND. DEN0M2 .NE. 0. .AND. 
1DEN0M3 .NE. 0. .AND. DENCM4 .NE. 0.) GO TO 30 
ZPC= I (WIN, JKWtN-5, J) 1/120. +( WIN-1, J ) +WIN-4, J ) )/24. ♦IWIN-2.J) 
1 *W I N-3 , J ) ) / 12 . ) 

RETURN 

30 CONTINUE 

DO 20 1=1,15 

20 Y I I ) =W I I , J) 

ZPC REPRESENTS THE 5TH DERIVATIVE OF THE FUNCTION 
THE USER MAY OEFINF ZPC ANALYTICALLY OR USE THE STATEMENT 
FUNCTION DERIV WHICH FINDS ZPC VIA THE 4 TH OIVICEC DIFFERENCE 
OF THE FUNCTION'S OER I V AT I VE , FOOT 
CONTINUE STATEMENT REPRESENTING THE EON NUMBER 
GO TO (1, 2, 3, 4, 5, 6, 7, 8, 9,10), J 

1 CONTINUE 

ZPC=DERI V (FN,FNI,FN2,FN3,FN4,DEN0M,DEN0M1 » DEN0M2, DEN0M3, DEN0M4) 
RETURN 

2 CONTINUE 

ZPC=DERI V IFN f FNl,FN2,FN3,FN4, CENOM.DE NON 1, DEN0M2.DEN0M 3.DENOM4) 
RETURN 

3 CONTINUE 

ZPC=OER IV (FN,FN1,FN2,FN3,FN4, DFNOM.DENOMl ,DEN0M2,DEN0M3,DEN0M4) 
RETURN 

4 CONTINUE 

ZPC=DER I V I FN.FNl ,F N2 ,FN 3 ,FN4 , DE NCM ,DENCM 1 .DEN0M2 , CEN0M3 , DENCM4 ) 
RETURN 

5 CONTINUE 
RETURN 

t CONTINUE 

RETURN 

7 CONTINUE 

RETURN 

6 CONTINUE 
RETURN 

9 CONTINUE 
RETURN 

10 CONTINUE 
RETURN 
END 
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SUBROUTINE RKUTTA 
IMPLICIT RE AL*8 (A-H.O-Z) 

COMMON/ZILCH/ FDOT ( 10 , 1 5 i , W ( 10 , 1 5 ) ,CST ( 15 ) , YYI 10,15),Z(15),T<10), 

1 ZIPI 15) .FOOTS! 10,15), PCER ,DELT , DELMAX , DEL .SAVE, SI , SAVE 1 , SII, 

2 I SKI PI 15) ,K,N, J, JMAX, I STEP, IDUBLE, 1CHECK.NQDBLE ,NONLIN,I TERAT 
01 MENSI ON A ( 4 ) 

00 10 1=1,4 
CALL FUNCT 

A ( I ) = OELT * FOOT ( K , J ) 

GO TO ( 1,1,3,10) , I 

1 I F l J . EQ. 1 ) T (K + l ) = T ( 1 ) +OFLT *.5 

2 W ( K+l , J ) = W(l,J) + A(I)*. 5 
K=K+ 1 

GO TO 10 

3 I F ( J .EO.l) T (K + l ) = T ( 1 ) +DELT 
VI ( K+ 1 , J ) = W( 1, J) + A ( I ) 

K=K+ 1 

10 CONTINUE 

N= I SK I P ( J ) +4 

W(N,J)=W(l,J)+(A(l)+2,*(A(2)+A(3))+A(4))/6. 

I F ( J .EO.l) T ( N ) = Til) +OELT 
H(K,J)=W(N,J) 

T(K)=T(N) 

CALL FUNCT 
FOOT(N,J)=FOOT(K,J) 

K= 1 

IF(ISKIPIl) .GE. 2) GO TO 50 
IF ( ICHECK .GE. 1) GO TO 15 

CALL STPCHK(W,FDOT,T,OELT,DEL, J, JMAX, I CHECK, N) 

SAVE 1 = 1 . 

RETURN 

15 I F ( ICHECK .EQ. 2) GO TO 20 

IF ( ICHECK .EQ. 1) CALL STPCK 1 ( W, FOOT ,T ,N, J, JMAX , I CHECK ) 

SAVE 1=1. 

RETURN 

20 I F ( J .LT. JMAX) RETURN 

CALL STPCK2(W,FD0T,T,SAVE1,N,J, JMAX, ICHECK) 

C IF ICHECK=3 DELT WAS TOO LARGE ,CALC WILL BE REDONE WITH CELT/2 
IF (ICHECK .NE. 3) GO TO 30 
ICHECK=0 
RETURN 

50 I F ( J .LT. JMAX) RETURN 

30 DC 25 J= 1 » JMAX 

W ( 1 , J ) =W ( N, J ) 

FOOT ( 1 , J ) =FOOT ( N, J ) 

ISKIP(J) = I SKI P { J )♦ l 

THE FOLLOWING STATMENT LOADS THE 9TH ARRAY LOCATION OF FD0T(9,J) 
AND W ( 9 , J ) SO THAT WHEN THE PRED-CORR TAKES OVER THE 9TH 
LOCATION WILL CONTAIN MEANINGFUL VALUES 
IF ( N .EQ. 8) W(N+L,J)=W(N,J) 

IF (N .EQ. 8) FOOTCN+l, J)=FDOT(N, J) 

100 FORMAT! 1H0,3D15. 8, 215) 

25 CONTINUE 

K= 1 

T(1)=T(N) 

J= JMAX 
RETURN 
END 
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SURROUTI ME STPCHK ( U , FOOT ,T ,OELT , DEL, J, JMAX , I CHECK* N ) 

IMPLICIT REAL*8 (A-H.O-Z) 

DIMENSION A ( 15) ,0 I 15 ) , W I 10 , 1 5 ) , FOOT ( 10, 15 ) . T { 10>,D(15),E(15) 
A(J)=W(N,J> 

8 ( J ) =FDOT I N » J ) 

D( J)=W( l, J) 

E I J ) =FOOT ( 1 » J ) 

I F ( J .LT. JMAX ) RETURN 
C=T ( N ) 

ICHECK= 1 
F=T ( 1 ) 

DELT=DELT*.5 

RETURN 

ENTRY STPCKK W, FOOT, T,N, J, JMAX, I CHECK) 

W ( 1 » J ) =W I N» J ) 

FDCTI 1, J)=FDOT(N, J) 

I F ( J .LT. JMAX) RETURN 
ICHECK=2 
T 1 1 ) = T ( N ) 

RETURN 

ENTRY STPCK2 (W, FOOT, T. SAVE l.N, J, JMAX, I CHECK) 

DO 1 1=1, JMAX 
CJHECK=DA0S(AI I) -WIN, I ) 1 

I F (CHECK .GT. 1.0-05*DABS (WIN, I ) ) ) GO TO 5 

1 CCNTINUE 
DELT=OEL 
ICHECK=0 
SAVE 1 = 0. 

RETURN 

5 CCNTINUE 

DO 21=1, JMAX 
W( 1, I )=D( I ) 

FOOT ( 1, l ) =E ( I ) 

2 CCNTINUE 
T ( 1 ) = F 
DEL=DELT 
ICHECK=3 
SAVE 1 = 1 . 

RETURN 

ENO 
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SUBROUTINE INTGRT 
IMPLICIT REAL*8 (A-H,0-Z) 

EXP!QQQ)=DFXP( OOQ ) 

ALOG!OQQ)=OLOG(QQO) 

ABS!QQQ)=DABS!QQQ) 

COMMON/ZILCH/ FOOT ! I 0 , 1 5 ) , W ( 10, 1 5 ) , CST I 15 ) , YY| 10, 15 1 , Z C 15 ) , TI 10) , 

1 ZIP( 15) , FOOTS! 10,15) ,PCER ,0E l T ,CELMAX ,DEL , S AVE , S I , S AVE1 , SII, 

2 I SKIP! 15) ,K,N, J, JMAX.I STEP, IDUBLE, ICHECK,NODBLE,NONL IN, ITERAT 
REAL MOO I FR 

DIMENSION PRFD! 15),C0RR( 15) 

22223 CONTINUE 

IF! I SK I P ! J ) .GT. A) GO TO 20 

C SAVE ORIGINAL PRED-CORR ERROR CRITERIA IN CASE OF LATER CHANGE 
PCERS=PCER 

C SET ITER =0 INITIALLY 

ITER=0 

160 CALL RKUTTA 
RETURN 
20 N=8 

K = N 

C FIFTH DERIVATIVE FOR THE ERROR TERM 

NZPC=N 
TTEMP=T ( N ) 

CALL DERIV5!TTEMP,ZPC,W,FD0T,T,J,NZPC) 

C PREDICTOR ! J ) IE FOR THE JTH EON 

PRED! J)=W!N,J)M0ELT/24. )*! 55.*FDOT ! N, J )-59. *FDOT! N- 1 , J) *37. * 
l FDOTIN-2 , J) -9. *FDGT! N-3 , J) ) ♦ ! 2 5 1 . / 720 . ) *DELT**5*ZPC 
MOD I FR=PRED I J ) 

W!N+1,J)=M001FR 

I F ( J . EO. 1 ) T!N*1)= TIN) +DFLT 

TTEMP=T!N+1) 

I TERAT=0 
600 CONTINUE 
N=S 
K = N 

CALL FUNCT 
NZPC=N 

CALL DERIV5I TTEMP,ZPC, W ,FDOT ,T , J , NZ PC ) 

N=8 

C CORRECTOR IJ) IE FOR THE JTH FQN 

CORR! J)= W(N,J)+(0ELT/2A. ) * ( 9.* FOOT ( N* 1 , J )•*• 19. *FDOT ( N, J )-5.*FD0T I N 
1-1, J )+FDOT !N-2, J » ) -! 19. /720. ) *DELT**5*ZPC 

C ITERATE ON CORRECTOR TO FIND CORRECT VALUE OF FUNCTION 

WTEMP=W!N+1, J) 

W(N+1,J)=C0RR(J) 

ITERAT=I TERAT+1 

IF! ITERAT .LT. 1000) GO TO 700 
IFIITER.EQ. 5) CALL DUMP 
ITER=ITER+1 

WRITE! 6,11000 PRED ( J ) ,CGRR ( J ) , FOOT ( 9, J ) , 

1TTEMP, J 
700 CONTINUE 

CHECK=AB S( W TEMP-W ( N+ 1 , J ) ) 

IF! CHECK .GT.1.E-07*ABS!W(N+1, J) ) )G0 TO 600 
11000 F0RMAT!1H0,AD16.8,I5> 

ZIP(J)= PRFD (J) -CORR(J) 

IF (CORR! J) .NE. 0.)GC TO 300 
Z( J)=0. 

GO TO 305 

300 Z!J)=PRED!J) /CORR ! J ) 
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305 CONTINUE 

CALL ERROR 

I F I J .LT. JMAX .AND. NONLIN .LE. 1) RETURN 
C THE IFS ON NONLIN SET UP THE ITERATION OF A SYSTEM OF NONLINEAR EGN 
IFINONLIN .EO. I) CALL SYSIT 
135 CONTINUE 

IFINONLIN .GE. 21 CALL SYS I T 1 1 6600, 0125, £ 135 ) 

125 CONTINUE 

C STEP SIZER CHECKS ERROR TERMS HERE 

C ERROR CRITERIA DETERMINES WHETHER TO DOUBLE OR HALF 

C SI AND SAVE CHECK STABILITY AND TRUNCATION ERROR RESP 
IF ( SI .GT. SII) GO TO 211 

C PCER IS THE ERROR CRITERIA DEPENDING ON THE PRED AND CORR. 

C PCER IS MADE STRONGER IF SI=0. I DFDY=0) SINCE IT IS THEN THE ONLY 

C ERROR CRITERIA FOR CHECKING 

IFISI .EQ. 0. .ANO. PCER .GT. .0001) PCER=.0001 
I F ( ABS ( SAVE I .LE. PCER) GO TO 50 
211 CALL STPSIZ 
IDUBLE=1 
NGDBLE=0 
GO TO 106 

50 IF (NODBLE .EO. 1) GO TO 105 
I F ( I STEP .LE. 5)GC TO 1C5 

IFISI .GT..5*SII.AN0. SAVE .GT. .5*PCER) GO TO 105 
CALL DOUBLE 

106 I F ( DEL .EO. DFLT ) GO TO 105 
DELT=DEL 
SAVE1=1. 

GO TO 110 
105 CONTINUE 

DO 115 J-l , JMAX 



DO 115 N= 1,8 

WIN, J)=WIN+1, 

J) 


I F ( J .EQ. 1 ) 

T( N ) =T I N*- 1 ) 


FDCTIN,J)=F00T(N+1,J) 

115 

CONTINUE 


116 

ISTEP=ISTEP*1 
J = JMA X 
CONTINUE 



SAVE 1=0. 

I 0U8LE= l 

IF l I STEP .GT. 

100) I S TE P = 5 


C RESET PCER IN CASE OF CHANGE 
PCER=PCERS 

C RFSET ITER IF CHANGE HAS OCCURREO 

IF 1 1 TER .GT. 0) I TER=0 

C THE FOLLOWING STATEMENTS CHECK THE STEP SIZE AND LIMIT DELT WHEN 
C DELMAX .GT. 0,IE., STEP LIMITER 
IF (DELMAX. EO. 0) GO TO 110 

I F I DELT .EO. DELMAX .OR. 2.*DELT .GT. DELMAX) N0DBLE=1 
110 CONTINUE 
22221 CONTINUE 
RETURN 
END 
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SUBROUTINE SYSIT 
IMPLICIT REAL*8 (A-H.O-Z) 

COMMON/ ZILCH/ FOOT! 10, 15 J ,W I 10, 15 ) ,CST ( 15 ) , YY( 10* 151 , Z ( 15 ) , T ( 10 ) , 

1 Z IP( 15 I , FOOTS! 10 ,15 ) ,PChR , DELT ,CELM AX, DEL .SAVE, SI , SAVE 1, S 1 1 . 

2 ISKIPI 15) ,K,N, J, JMAX, ISJEP, IOUBLE, ICHECK,NODBLE,NONLIN, I TER AT 
DIMENSION YDOTI 15) 

C SYSIT PERFORMS ITERATION ON SYSTEMS OF NONLINEAR DIFFERENTIAL EQN 

1 IFINONLIN .EO. 3) NONLINM 

2 DO 3 J= l , JMAX 

YDOT( J)=FD0T(9,J) 

3 CALL FUNCT 

DO 5 1=1, JMAX 

CHECK=OABS I YDOT ( I )-FDOT (9,1 ) ) 

IF! CHECK .GT. I . E-07*DABS I FOOT I R , 1)) ) GC TO 1 

5 CONTINUE 

IFINONLIN .EC. 3) GO TO 6 
J=C 

IFINONLIN .60. 4) RETURN 3 

NCNL I N=2 

RETURN 

ENTRY SYSIT1I*,*,*) 

J = J+1 

IFIJ .LE. JMAX) RETURN l 

C ITERATION TO CHECK COMPATABILITY OF FOOT WITH NEW VALUFS FOR W 
NCNL I N=3 
GC TO 2 

6 CONTINUE 
NCNL I N= 1 
J= JMAX 
RETURN 2 
END 
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c 

c 


25 

20 

100 


11 


2 


3 


4 

1 

30 


SUBROUTINE STPSIZ 

IMPLICIT RE AL*8 IA-H.O-Z) Tlin . 

COMMON/ Z I LCH/ FOOT ( 10, 1 5 ) * W 1 10, 1 5 » , CST ( 1 5 1 , YY 1 10, 15 1 * Z 115 j » T ( 10 * 

1 ZIP! 15) , FOOTS! 10,15) ,PCER ,0E L T .CELMAX ,DEL , S AVE , S I , S AVE1 , S 1 1 . 

2 I SKIP! 15) ,K,N, J, JKAX.ISTEP, 1 DUPLE, I CHECK , NODBLE , NONL 1 N, ITERAT 
DIMENSION WS!8)*TX<8)*FS(8)*A!1Q)*B(10) 

DEL=DELT* . 5 


I S = 2 

10RDER=4 

N0PRNT=0 

IF OOUBLING AND HALVING OCCURS ON SUCCESSIVE STEPS USE ACTUAL 

VALUES FOR W!7,J) AND FD0T!7,J) 

IF t IDU8LE .EO. 1) GU TO 20 
DO 25 J= 1 , J MA X 
W<4, J)=YY!4, J) 

W< 5 , J ) =YY ! 5 , J ) 

W ( 6 , J ) =YY ! 6 , J ) 

W ( 7 , J ) = YY ( 7, J ) 

FDCT!2,J)=FD0TS!2,J) 

FDOT(3,J)=FDOTS!3,J) 

FOCT!A,J)=FDOTS(A,J) 

FOOT ( 5,J )=FDOTS! 5,J) 

FOCT!6,J)=FDOTS!6,J) 

FOOT ( 7, J ) =FDOTS ! 7 , J ) 

FDCT19,J)=F00TS!9,J) 

CONTINUE 
GO TO 30 
CONTINUE 
DO 1 J= 1 , JMAX 
FORMAT! IX, 8015. 8) 

DO 11 1=5,8 
WS! I - A ) = W ( I , J) 

FS ( I-A ) =FOOT ( I, J) 


(XI,WS,TX,IMAX,IS,L,ICRDER,UX,NOPRNT) 

! XI ,FS,TX, IMAX, IS,L, IORDER,UX,NOPRNT) 


W(M,J)=W!MM,J) 

FDCT(M,J)=FDOT!MM,J) 

CONTINUE 

MM=1 

DO A M=3 ,7,2 
W ! M » J ) = A 1 MM ) 

FOOT ! M , J ) =8 ! MM) 

MM=MM+1 

CONTINUE 

CONTINUE 

CONTINUE 

I0RDER=2 

I MAX= 8 

DO 5 K=5 , 7 


TX! 1-4) =T< I ) 
CONTINUE 
DC 2 L= 1 , 3 
Xl=TX ! L ) +DEL 
CALL SI N TP 
AIL) =UX 
CALL SINTP 
BID =UX 
CONTINUE 
MM = 4 

DC 3 M=2,6,2 
mm=mm+i 


D-29 



X I=T ( K ) +DEL 

CALL SINTP ( X I » T , T , IMAX, l S, K, 10RDER ,UX ,NOPRNT ) 
TX(K)=UX 

5 CONTINUE 
MM = 5 

1=4 

DO 6 K= 1 « 7 

GC TO (15, 16, 15, 16,15, 16, 15), K 

15 T(K)=TX(I) 

1 = 1 + 1 

GO TO 6 

16 T(K)=T(MM) 

MM=MM+1 

6 CONTINUE 
J= JMAX 
RETURN 
ENTRY DOUBLE 

C STEP SIZE =DEL T*2 
DEL=OELT*2. 

DO 55 J= 1 , JMAX 

I DU8LE=2 

YY ( 4, J ) =VI ( 4, J ) 

YY(5,J)=W(5,J) 

YYI 6 , J ) =W ( 6 , J ) 

Y Y ( 7 , J ) =W ( 7, J ) 

FOOTS(2,J)=FDOT(2,J) 

FD0TS(3,J) =FDOT ( 3 , J ) 

FDOTS(4,J)=FOOT(4,J) 

FOOTS ( 5 , J ) =FDOT ( 5 , J I 
FDOTS( 6,J) =FDOT( 6 , J ) 

FD0TS(7,J) =FDOT 1 7 , J ) 

FDCTS(9, J J=F00T(9, J) 

DO 55 1=1,3 
GO TO (65,70,75), I 
65 N=7 

NN = 6 

GO TO 80 
70 N=6 

NN=4 

GO TO 80 
75 N=5 

NN=2 

80 CONTINUE 

W( N, J )=W ( NN, J ) 

I F ( J . EQ. 1 ) T(N)= T ( NN ) 

FDOT(N, J)=FOOT(NN, J) 

55 CONTINUE 

I STEP=i 
J= JMAX 
RETURN 
ENC 
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SUBROUTINE SINT^ (A,B,C , I M AX , I S , l , I OROER, D *NOPRNT ) 

IMPLICIT REAL*8(A-H,0-Z) 

DIMENSION B( ICO) ,C( ICO) 

DIMENSION X!100) ,YI 100),NUMI 20),DEN0M< 20, 20),U(20) 

REAL NUM.NUMT 
X I = A 

DO 60 K=l , IMAX 
Y I K ) =B { K ) 

X I K ) =C ( K ) 

60 CCNTINUE 

I F ( I OROER .LE. 1) I0RDER=2 
SEARCH AND INTERPOLATE ROUTINE 
IS=0 FOR BOTH SEARCH AND INTERP 
l S= 1 FOR SEARCH ONLY 
l S = 2 EOR INTERP CNLY 

N0PRNT=3 SUPRESSES THE WRITTING OF THE 

DI AGON I ST I C MESSAGE , NOPRNT RETURNS A VALUE OF -1,0, OR L DEPENDING 
ON WHETHER THE INTERPOLATION TOOK PLACE OUTSIDE THE TABLE , LOWER 
(-1 ) ,UPPERENO( 1) ,0R INSIDE THE TABLE ( 0 ) 

IF ( I S .FO. 2 ) GO TO 20 
1 = 1 

I F ( X I .LT. XI IMAX) .AND. XI .6T. XII)) GO TO 10 
IFIXl .GT. X( I) ) GC TO 15 
1 = 1 

GO TO 20 
15 I = I MAX 

GO TO 20 

10 IFIXl .LE. XU + 1) .AND. XI .GT.X(D) GO TO 20 
1 = 1 + 1 
GO TO 10 

20 CCNTINUE 

IFIIS .EG. 1) RETURN 

C TABLE WITHIN A TABLE INDEX LOCATION 
I T= 1 

IFIIS .EG. 2 .AND. I .GT.1)IT=I 
I F ( IOROER .GT. IMAX) I ORDER= IMAX 
I OROR S= I ORDER 

22 I F ( ( l + I ORDER- 1 ) .LE. IMAX) GO TO 21 

I0RDER=I0RDER-1 
I F ( I ORDER .NE. 1) GO TO 22 
A 6 I F ( NOPRNT .EG. 3) GC TO A5 

WRITE (6, 101) Y( IT) ,Xf IT) , Y ( I MAX ) , X ( I MAX ) , X I , I MAX 
101 FORMAT! 1H0, 'AN ERRONEOUS VALUE HAS BEEN SUPPLIED TO SINTRP, IE XI 

1 IS OUTS IDF THE TABLE BEING USED , Y1 = ' , C 15 . 8 , * X 1 = * , C 1 5 . 8 , / , • YMAX = 

2 • ,015. 8, ' XMAX= • , 015.8,’ XI=* ,D15. 8,' IMAX = ' ,1 5, /, • THE APPROPRIATE 
3END POINT IS USED FOR THE INTERPOLATED VALUE') 

A 5 CCNTINUE 
I I = I MAX 

I F ( X I .LT. X ( I M AX ) ) I I = 1 
UX = Y ( I I ) 

IQRDER= IORDRS 
NOPRNT=- 1 

IF! II .FO. IMAX) N0PRNT=1 

D = UX 

RETURN 

21 IFIXl .LT. X! IT ) ) GO TO A6 

C SAVE THE ORIGINAL INDEX I TO RETURN TO CALLING PROGRAM 

I RETRN= I 

C THE INDEX I REPRESENTS THAT POINT IN THE TABLE WHERE THE 
C POLYNOMIAL WAS STARTED 
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I F ( I ORDER . NF . IORDRS) I = I-t IORDRS- 1 ORDER I 
NTH ORDER INTERPOLATOR 
0 CONTINUE 

RETURN IORDER TO ORIGINAL VALUE 
I ORiJER= I ORDRS 
NOPRNT = 0 
l SAVE = 1 

DO 25 1 1 = 1 1 1 ORDER 

I = I SAVE* II — 1 
NUMI I I ) = XI - X C I) 

DO 25 LL=1 » IORDER 
L=LL + I SAVE-1 

OENOMI 1 1 » LL ) =X( I)-X(L) 

SET DIAGONAL =1 SO THAT DIVISION BY 0 DOES NOT OCCUR 
I F ( I .EO. L) DENOMI II ,LL)=1. 

25 CONTINUE 

00 26 1=1, IORDER 
TERM =1. 

DO 27 L=1 ,IOROER 

I I =L 

IF ( I .NE. L) GO TO 30 
NUNT = NUM C I I ) 

NUN ( I I ) = 1 . 

3 0 TERM=TERM*l NUMI I I ) /OENOMI I »L I ) 

IFII .EQ. L) NUMI 1 1 ) =NUMT 
27 CONTINUE 

U I I ) = TERM 

26 CONTINUE 
1 = l SAVE 

C VALUE OF INDEP VAR AT XI 
UX = 0. 

00 AO K= 1, IORDER 
UX=UX+UIK)*Y< I) 

1 = 1 + 1 

AO CONTINUE 

I=IRETRN 
D = L X 
RETURN 
END 

.. DATA DMS.N.GCBIR 
/* 


NASA— GSFC 
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